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The  use  of  fiber  reinforced  composite  laminates  in  engineering  applications  has  been 
increasing  rapidly.  Along  with  this  increase  has  come  a  rapid  development  in  the  analysis 
techniques  to  accurately  model  internal,  as  well  as  gross  plate  behaviors.  Many  improve¬ 
ments  to  laminated  plate  theory  have  been  developed  in  the  push  for  better  analysis 
techniques.  Improvements  began  with  the  application  of  Mindlin-Reissner  shear  deforma¬ 
tion  theory  followed  by  higher  order  theories  and  discrete  layer  theories.  With  the  drive  for 
more  accurate  modeling,  the  cost  has  been  increased  complexity  and  computational  time. 
Some  of  the  higher  order  techniques  lend  themselves  well  to  simplification,  but  in  doing  so 
they  complicate  the  finite  element  analysis  by  creating  a  C1  continuity  requirement.  The 
purpose  of  this  work  is  to  provide  accurate,  yet  computationally  efficient,  improvements 
to  the  analysis  of  composite  laminates. 

One  portion  of  this  work  shows  that  the  higher  order  extensions  to  the  first  order 
shear  deformation  theory  still  do  not  correctly  model  the  physics  of  the  laminated  plate 
problem.  Results  show  that  the  first  order  theory  cam  provide  as  good,  if  not  better,  re¬ 
sults  with  the  proper  sheax  correction  factor.  This  work  uniquely  implements  a  Predictor 


Corrector  technique  into  the  finite  element  method  to  accurately  calculate  the  shear  cor¬ 
rection  factors.  The  technique  provides  excellent  results  with  a  simple  Mindlin  type  plate 
element. 

The  second  part  of  this  research  develops  two  new  finite  elements  which  approximate 
C  continuity  through  the  use  of  a  least  squares  technique.  These  Least  Squares  elements 
can  be  used  to  take  advantage  of  the  displacement  field  simplification  techniques  which,  up 
until  now,  have  seriously  complicated  the  finite  element  application.  The  implementation 
of  the  elements  are  demonstrated  using  a  piecewise,  simplified  third  order  displacement 
field.  The  Least  Squares  elements  should  prove  to  be  useful  tools  in  any  finite  element 
application  where  C1  continuity  is  required. 

The  final  portion  of  this  work  presents  a  study  into  the  effects  of  stacking  sequence, 
boundary  conditions,  pre-stress  and  plate  aspect  ratios  on  the  fundamental  frequency  and 
buckling  loads  of  laminated  plates. 
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SUMMARY 

The  use  of  fiber  reinforced  composite  laminates  in  engineering  applications  has  been 
increasing  rapidly.  Along  with  this  increase  has  come  a  rapid  development  in  the  analysis 
techniques  to  accurately  model  interned,  as  well  as  gross  plate  behaviors.  Many  improve¬ 
ments  to  laminated  plate  theory  have  been  developed  in  the  push  for  better  analysis 
techniques.  Improvements  began  with  the  application  of  Mindlin-Reissner  shear  deforma¬ 
tion  theory  followed  by  higher  order  theories  and  discrete  layer  theories.  With  the  drive  for 
more  accurate  modeling,  the  cost  has  been  increased  complexity  and  computational  time. 
Some  of  the  higher  order  techniques  lend  themselves  well  to  simplification,  but  in  doing  so 
they  complicate  the  finite  element  analysis  by  creating  a  C1  continuity  requirement.  The 
purpose  of  this  work  is  to  provide  accurate,  yet  computationally  efficient,  improvements 
to  the  analysis  of  composite  laminates. 

One  portion  of  this  work  shows  that  the  higher  order  extensions  to  the  first  order 
shear  deformation  theory  still  do  not  correctly  model  the  physics  of  the  laminated  plate 
problem.  Results  show  that  the  first  order  theory  can  provide  as  good,  if  not  better,  re¬ 
sults  with  the  proper  shear  correction  factor.  This  work  uniquely  implements  a  Predictor 
Corrector  technique  into  the  finite  element  method  to  accurately  calculate  the  shear  cor¬ 
rection  factors.  The  technique  provides  excellent  results  with  a  simple  Mindlin  type  plate 
element. 

The  second  part  of  this  research  develops  two  new  finite  elements  which  approximate 
C1  continuity  through  the  use  of  a  least  squares  technique.  These  Least  Squares  elements 
can  be  used  to  take  advantage  of  the  displacement  field  simplification  techniques  which,  up 
until  now,  have  seriously  complicated  the  finite  element  application.  The  implementation 
of  the  elements  are  demonstrated  using  a  piecewise,  simplified  third  order  displacement 
field.  The  Least  Squares  elements  should  prove  to  be  useful  tools  in  any  finite  element 
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application  where  C1  continuity  is  required. 

The  final  portion  of  this  work  presents  a  study  into  the  effects  of  stacking  sequence, 
boundary  conditions,  pre-stress  and  plate  aspect  ratios  on  the  fundamental  frequency  and 
buckling  loads  of  laminated  plates. 


CHAPTER  I 


A  REVIEW  OF  THE  LITERATURE  AND 
DEVELOPMENTS  IN  LAMINATED  PLATE 

THEORIES 


1.1  Introduction 

Laminated  fiber  reinforced  composite  materials  have  provided  engineers  with  the  ability 
to  design  and  build  structures  as  never  before.  The  use  of  composites  has  been  growing 
rapidly  over  the  past  twenty  years  and  is  continuing  to  do  so  at  an  increased  rate.  Early  in 
their  existence,  their  use  was  primarily  associated  with  spacecraft  and  aircraft  because  of 
their  high  strength  to  weight  ratios,  in  spite  of  their  high  cost.  Recently  however,  reduced 
manufacturing  costs  are  making  composites  attractive  to  many  other  industries.  Compos¬ 
ites  are  now  being  used  for  automobiles,  sporting  goods,  pressure  vessels  and  a  multitude 
of  other  applications.  Composite  materials  will  eventually  be  able  to  benefit  virtually  any 
engineering  application  because  of  their  design  advantages.  Today’s  technology  has  only 
begun  to  realize  the  resource  that  is  becoming  available  in  the  composite  material  world. 
The  engineer  has  the  ability  to  not  only  design  directional  strength,  but  also  thermal  and 
electrical  conductivity,  radar  absorption,  thermal  expansion,  fracture  characteristics  and 
stiffness,  to  only  mention  a  few  parameters.  As  research  into  composite  materials  con¬ 
tinues,  more  and  more  of  these  design  parameters  will  be  developed,  and  more  and  more 
applications  will  arise.  In  reality,  the  composite  material  science  is  probably  in  its  very 
infancy,  and  as  it  continues  to  grow,  so  must  the  ability  to  perform  accurate  engineering 
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analyses. 

The  engineering  analysis  of  composite  materials  is  in  itself  a  relatively  new  field  and  has 
just  begun  to  grow.  The  mathematical  modeling  of  the  mechanics  of  composite  materials 
dates  back  only  thirty  years  ago  when  classical  laminated  plate  theory,  as  we  know  it 
today,  was  developed  by  Reissner  and  Stavsky  (1961)  [121].  It  remains  today  as  the  main 
tool  available  to  the  practicing  engineer.  However,  as  the  field  grows  so  will  its  complexity, 
and  classical  laminated  plate  theory  (CLPT)  will  not  be  a  sufficient  analysis  tool.  As  the 
field  grows,  more  accurate  and  efficient  modeling  techniques  must  be  developed.  The 
inherent  complex  nature  of  composite  laminates  often  necessitates  complex  mathematical 
models.  Unfortunately,  complex  models  are  difficult  to  implement  in  practical  engineering 
analysis,  so  the  need  for  accurate,  yet  efficient,  methods  will  remain  high.  No  matter  how 
accurate  or  simple  a  mathematical  model  is,  it  has  very  limited  engineering  applicability  if 
it  cannot  be  applied  to  general  shapes  and  boundary  conditions.  The  finite  element  method 
is  the  tool  which  is  generally  used  to  achieve  this  capability.  However,  the  finite  element 
implementation  of  new  mathematical  formulations  can  be  difficult  and  the  end  product 
is  not  always  useful.  Accuracy  in  the  finite  element  method  many  times  corresponds  to 
increased  computational  costs.  For  example,  a  recent  article  by  Jing  and  Liao  (1989) 
[39]  proposes  a  new  element  which  gives  excellent  results  for  laminated  composites.  The 
element  is  employed  in  each  layer  of  the  laminate.  Thus,  each  layer  is  modeled  by  a  twenty- 
node  mixed  field  hexahedron  with  three  degrees  of  freedom  at  each  node  and  fourteen  stress 
parameters.  One  can  see  that  for  a  laminate  with  a  moderate  number  of  layers  the  analysis 
can  quickly  become  numerically  intractable. 

Based  upon  the  above  discussion,  we  see  that  the  future  calls  for  not  only  increased 
understanding  and  more  complex  mathematical  modeling  of  composite  materials,  but  also 
for  fresh  ideas  and  approaches  on  how  to  effectively  and  economically  model  laminate 
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behavior.  It  is  hoped  that  this  work  will  present  some  novel  approaches  in  the  analysis 
methods  of  composite  materials  which  will  provide  simple,  yet  powerful,  tools  to  be  used 
in  engineering  design  analysis.  In  addition,  it  may  possibly  initiate  a  new  methodolog* 
for  future  work  in  laminated  composite  plates  and  shells. 

1.2  A  Brief  Review  of  Basic  Plate  Theory 

The  mathematical  analysis  of  plates  has  been  a  much  studied  area  in  the  engineering  world 
for  many  years.  The  use  of  plates  as  major  structural  components  has  driven  researchers 
to  find  a  way  to  accurately  predict  their  behavior  from  a  static,  dynamic  and  stability 
point  of  view.  The  first  major  achievements  in  modem  engineering  plate  analysis,  as 
stated  by  McFarland  et  al  (1972)  [74],  were  begun  in  the  early  1800's  and  are  accredited 
to  Cauchy,  Poisson,  Navier,  Lagrange  and  Kirchhoff.  However,  the  development  of  what 
we  know  today  as  classical  plate  theory  (CPT)  is  generally  attributed  to  Kirchhoff  [55]  for 
his  work  in  1850. 

In  CPT  certain  assumptions  are  made  simplifying  the  problem  to  one  that  is  more 
easily  solved.  The  Kirchhoff  assumptions,  as  they  are  sometimes  called,  parallel  the  ideas 
behind  simple  beam  theory.  We  first  assume  that  a  normal  to  the  midplane  of  the  plate 
before  deformation  remains  normal  and  inextensionable  after  deformation.  Also,  we  as¬ 
sume  that  normal  stresses  in  the  transverse  direction  to  the  plate  are  small  compared 
with  the  other  stresses  and  can  be  neglected.  The  geometry  of  the  deformation  is  shown 
in  Figure  1.1.  One  can  see  that  the  in-plane  displacements  are  composed  of  a  translation 
and  a  rotation.  They  can  be  written  as: 

dw„ 

u  =  ti0  z  ~T 
ox 

dw0 
dy 


v 


V0  ~  Z 


(1.1) 


4 


Z 


Figure  1.1:  Deformation  Geometry  for  CPT 


w  =  u>0 

The  Kirchhoff  assumptions  are  valid  for  many  cases,  and  accurate  results  can  be 
achieved  with  them  for  engineering  problems.  Problems  are  restricted  to  thin  plates  free 
from  any  large  transverse  loads.  However,  there  is  an  important  concept  to  remember 
when  working  with  the  Kirchhoff  assumptions.  One  must  remember  that  in  assuming 
that  the  normals  to  the  midplane  remain  normal  after  deformation,  one  does  not  preclude 
transverse  stresses  '.  Just  as  in  beam  theory,  it  means  that  the  additional  deformation 
caused  by  these  stresses  is  negligible.  This  is  a  valid  assumption  as  long  as  the  shear 
rigidity  for  the  transverse  strain  is  on  the  same  order  of  magnitude  as  the  elastic  modulus, 
which  is  the  case  for  most  isotropic  engineering  materials. 

The  next  major  advancement  in  plate  theory  was  the  logical  step  to  include  the  effects 
of  transverse  shear  deformation  into  the  governing  equations.  Including  transverse  shear 

allows  the  normals  to  the  midplane  to  deform.  The  work  in  this  area  closely  parallels 

’Throughout  this  work,  'transverse  stresses’  and  ‘transverse  strains’  will  imply  the  shear  components 
only  ,  and  not  the  normal  components  (unless  otherwise  specified). 
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that  done  in  beam  theory.  Transverse  shear  deformation  effects  are  included  in  going  from 
Bernoulli- Euler  beam  theory  to  Timoshenko  beam  theory.  The  inclusion  of  transverse 
shear  into  plate  theory  has  taken  many  forms  and  was  proposed  in  many  different  ways  by 
several  investigators.  In  a  survey  article  Reddy  (1985)  [114]  presented  a  brief  account  of 
the  development  in  this  area.  It  appears  that  work  to  include  transverse  shear  effects  into 
plate  theory  was  first  published  by  Basset  [8]  in  1890,  followed  by  Reissner  (1945)  [118], 
Hencky  (1947)  [30],  Hildebrand  (1949)  [32],  Reissner  (1947)  [119]  and  Mindlin  (1951)  [76]. 
Today  the  development  of  shear  deformable  plate  theory  (SDPT)  is  sometimes  categorized 
as  a  Reissner- Mindlin  plate  theory.  The  difference  being  that  Reissner  used  a  stress  based 
approach  and  Mindlin  used  a  displacement  based  approach,  as  did  Basset,  Hildebrand  and 
Hencky. 

In  SDPT  the  Kirchhoff  assumption  that  the  normals  to  the  midplane  remain  normal 
after  deformation  is  removed.  The  result  is  that  the  displacements  in  the  u  and  v  direc¬ 
tions  are  no  longer  constrained  to  be  equal  to  the  rotation  of  a  midplane  normal.  Such 
a  deformation  is  depicted  in  Figure  1.2.  In  its  general  form,  this  deformation  can  be 


6 


represented  as  a  power  series  expansion  in  z,  with  the  number  of  terms  carried  in  the 
expansion  being  determined  by  the  order  of  the  theory  desired.  In  the  initial  work  with 
Reissner-Mindlin  SDPT,  the  displacements  are  assumed  to  be  of  the  form 


U  =  U0  +  Z<f)x 

(1.2) 

V  =  Vo  +  Z<f)y 

(1-3) 

w~w0  (1.4) 

where  u0,  vQ,  w0,  4>x  and  <f>y  are  all  functions  of  x  and  y.  Here,  the  u  and  v  displacements 
are  linear  in  z,  while  w  is  constant  with  respect  to  z.  Thus,  the  deformed  normal  would 
maintain  a  straight  line  appearance  but  be  inclined  to  the  midplane.  For  obvious  reasons, 
the  SDPT  with  this  form  of  displacements  will  be  referred  to  as  the  first  order  theory.  Just 
as  a  point  of  comparison,  if  the  Kirchhoff  assumption  is  invoked  on  these  displacements 
we  have 


(1-5) 

(1.6) 


which  reduces  the  deformation  back  to  a  translation  and  the  same  rotation  as  in  eqns 


(  1.1).  This,  in  effect,  couples  4>x  and  4>y  to  w.  In  SDPT  wa  is  uncoupled  from  <f>x  and  <f>y 
creating  the  required  deformation  through  the  thickness.  Substitution  of  eqns  (  1.2)-(  1.4) 


into  the  linear  strain  displacement  equations  gives,  for  the  six  strains: 


(1.7) 


(1.8) 
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Cz  = 


dw 


=  0 


dz 
du  dw 


7 «  =  -3“  +  -3-  =  <£r  +  ty,r 

£72  £71 

dt7  dl£7 

'>'•*  =  fe  +  ^  =  *» +  ”•* 

du  dv  , 

Try  —  "t"  —  7^o,y  ~b  ^o,x  d"  z( <pr,y  "h  <Py,x)  —  To  “I”  ^^-xy 


(1.9) 

(1.10) 

(1.11) 

(1.12) 


(Note  the  use  of  commas  denoting  partial  differentiation.) 

The  results  in  eqns  (  1.7)-(  1.12)  show  an  important  aspect  of  the  first  order  theory. 
This  is  the  fact  that  the  assumed  displacement  fields  for  u  and  v  presuppose  constant  values 
of  transverse  shear  strains  in  the  thickness  direction.  Hence,  for  a  homogeneous  plate  the 
shear  stress  is  also  constant  through  the  thickness.  Again  parallel  to  beam  theory,  this 
obviously  cannot  be  the  case,  as  the  top  and  bottom  surfaces  should  be  free  of  any  surface 
tractions  (for  the  free  vibration  case).  The  result  of  this  violation  is  that  the  amount 
of  transverse  strain  energy  is  overestimated,  and  the  plate  model  is  actually  more  stiff 
than  it  should  be.  To  remedy  this,  a  shear  correction  factor  is  normally  introduced.  The 
correction  factor  multiplies  the  transverse  shear  stiffness  coefficients,  thus  reducing  the 
stiffness  in  these  directions.  For  an  isotropic  material  the  coefficients  could  be  introduced 


as 


axz  =  k-G-yxz  (1.13) 

ayz  =  k-Gjyz  (1.14) 

The  value  for  k  can  be  shown  to  be  dependent  upon  the  cross  section  of  the  beam  or 
plate.  For  an  isotropic  material  with  a  rectangular  cross  section  the  value  becomes  5/6. 
This  factor  is  due  to  the  fact  that  the  exact  solution  is  parabolic.  This  shortcoming  of  the 
first  order  theory  has  not  proven  to  be  of  any  detriment  to  the  results  obtainable  using 
it.  With  the  correct  shear  correction  factor,  the  theory  provides  very  accurate  results 
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for  the  gross  plate  behavior.  It  is  understandable  that  the  integral  average  of  transverse 
stresses  can  be  predicted  accurately,  but  their  distribution  through  the  thickness  cannot. 
Also,  the  assumed  displacements  result  in  in-plane  stresses  which  are  antisymmetric  about 
the  midplane  of  the  plate.  Thus,  inaccurate  results  may  be  achieved  for  cases  where  the 
loading  conditions  may  preclude  such  a  solution. 

1.3  Further  Developments  in  Shear  Deformable  Plate  The¬ 
ory 

Since  the  development  of  SDPT  many  researchers  have  published  works  on  variations  of 
the  Reissner-Mindlin  approach.  Since  1957  many  higher  order  displacement  fields  have 
been  used  in  an  attempt  to  achieve  more  accurate  results  and  to  eliminate  the  shortcom¬ 
ings  of  the  first  order  theory.  The  desire  has  mainly  been  to  eliminate  the  need  for  a  shear 
correction  factor.  The  landmark  article  into  improving  first  order  SDPT  was  published  by 
Lo,  Christensen  and  Wu  (1977)  [67].  In  their  article  the  authors  review  some  of  the  differ¬ 
ent  higher  order  displacement  fields  which  have  been  used  over  the  years.  Most  significant, 
however,  is  that  they  themselves  present  a  higher  order  theory  using  displacements  of  the 
form: 


u  =  uQ  +  z<t>x  4-  z2ipx  +  z3(x 

v  =  va  + z<j>y  +  z2ipy  +  z3(y  (1.15) 

w  =  w0  +  z<j>z  +  z2tpz 

In  their  work  Lo  et  al  show  that  this  form  of  a  displacement  field  alleviates  the  necessity 
for  a  shear  correction  factor  for  homogeneous  plates,  and  that  one  can  get  better  results  for 
cases  with  certain  loading  conditions  which  result  in  non-antisymmetric  in-plane  stresses. 
They  also  briefly  discuss  some  of  the  other  forms  of  higher  order  displacement  fields  which 
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have  been  studied.  The  end  result  appears  to  be  that  the  required  form  of  the  displacement 
field  is  closely  tied  to  several  items.  The  choice  of  the  displacement  field  depends  on  the 
type  of  problem  being  solved,  which  variables  are  needed  as  a  result  of  the  analysis  and 
the  required  level  of  accuracy  for  those  variables. 

For  gross  plate  behavior,  and  accurate  in-plane  stresses,  the  first  order  theory  has 
provided  very  satisfactory  results  over  the  past  forty  five  years.  At  the  cost  of  increased 
complexity,  the  higher  order  approaches  have  not  been  applied  in  any  extent  to  isotropic 
plate  theory.  In  fact,  it  was  not  until  the  advent  of  laminated  composite  plates  that  higher 
order  approaches  were  considered  to  any  extent  at  all.  In  the  above  mentioned  work  by 
Lo  et  al,  their  third  order  SDPT  was  developed  and  demonstrated  for  isotropic  plates,  but 
it  was  immediately  applied  to  laminated  plates  [68],  the  real  motivation  for  their  work. 
It  is  the  fact  that  first  order  SDPT  has  many  disadvantages  when  it  comes  to  laminated 
composite  plates  that  has  prompted  the  search  for  an  improvement  in  the  analysis  of 
such  structures.  This  search  has  lead  to  the  multitude  of  higher  order  applications  of 
SDPT  that  cam  be  found  in  the  literature  over  the  past  twenty  yeairs.  However,  before  we 
cam  accurately  amd  efficiently  apply  the  theories  available  to  us,  we  must  first  study  the 
fundamentals  of  the  problem. 

1.4  Laminated  Fiber  Reinforced  Composite  Plates 

Research  into  the  analysis  of  laminated  plates  begam  in  the  late  1950’s  by  Stavsky  (1959) 
[133]  amd  Lekhnitsky  (1968)  [62]  (originally  published  in  Russiam  in  1957),  but  classical 
laiminated  plate  theory  (CLPT)  as  we  know  it  today  is  credited  to  Reissner  and  Stavsky 
[121]  for  their  work  in  1961.  The  development  of  the  CLPT  equations  will  be  developed 
later  in  Section  2.2.1,  so  for  now  we  will  concentrate  on  understanding  the  physics  of  the 
deformation  of  laminated  fiber  reinforced  composite  plates. 
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1.4.1  The  Lamina 

The  geometry  of  a  lamina  is  shown  in  Figure  1.3.  The  material  properties  and  constitutive 
equations  of  the  lamina  are  described  originally  with  respect  to  the  material  coordinates, 
the  x\,X2  and  X3  axis,  and  then  transformed  into  the  global  coordinates  through  a  simple 
transformation.  The  details  of  this  will  be  presented  in  Section  2.2.1.  At  this  point  it  is 
important  to  understand  the  material  properties  of  the  lamina  in  Figure  1.3.  The  elastic 
modulus  is  generally  much  higher  in  the  x\  direction  than  in  the  12  direction.  Typical 
values  for  E\/E<2  can  be  on  the  order  of  40.  In  addition,  and  perhaps  more  importantly, 
the  values  for  shear  rigidity  are  small  compared  to  the  elastic  modulus.  A  typical  value 
for  G23/E2  is  0.35,  while  G\z/E\  can  be  less  than  0.01.  These  big  discrepancies  in  rigidi¬ 
ties  make  shear  deformation  considerations  very  important  in  the  analysis  of  composite 
materials.  Early  investigators  found  that  they  could  no  longer  neglect  the  contribution  of 
transverse  shear  deformation  to  the  overall  deformation,  even  in  relatively  thin  laminates. 
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Figure  1.4:  Laminate  Deformation  Geometry 

(See  Reissner  (1945)  [118],  Mindlin  (1951)  [76],  Whitney  (1969)  [141]  and  Whitney  and 
Pagano  (1970)  [144].)  These  properties  of  a  lamina  are  compounded  when  multiple  lamina 
are  stacked  to  form  a  laminate. 

1.4.2  The  Laminate 

The  real  advantage  of  using  composite  materials  in  the  design  of  structures  is  realized 
when  multiple  layers  are  stacked  with  varying  ply  angles.  This  allows  the  engineer  to 
tailor  the  material  properties  to  fit  a  specific  purpose.  This  process,  while  providing  great 
design  capabilities,  creates  great  difficulty  in  accurate  analysis  of  the  plate  parameters.  A 
laminated  plate  may  appear  to  behave  externally  like  an  isotropic  plate,  but  internally  it 
is  totally  different.  Figure  1.4  depicts  the  cross  section  of  a  laminated  plate. 

As  discussed  above  in  Section  1.4.1,  transverse  shear  effects  are  extremely  important 
in  the  analysis  of  composite  materials.  The  inclusion  of  transverse  shear  effects  in  the 
study  of  composite  plates  was  first  done  by  Yang,  Norris  and  Stavsky  (1966)  [145]  in  their 
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work  on  elastic  wave  propagation  in  heterogeneous  plates.  This  first  work  was  nothing 
more  than  the  application  of  a  Mindlin  displacement  form,  as  in  eqns  (  1.2)-(  1.4),  to  the 
classical  laminated  plate  theory  equations.  This  approach,  despite  its  simplicity,  has  been 
the  mainframe  for  the  analysis  of  transverse  shear  effects  in  fiber  reinforced  composite 
laminated  places.  In  fact,  commercial  analysis  codes  used  in  industry  rely  upon  this  first 
order  SDPT.  This  theory  has  been  used  in  one  form  or  another  extensively  over  the  past 
twenty  eight  years.  Most  of  the  advances  and  improvements  in  the  analysis  methods  since 
1968  have  been  extensions  in  one  form  or  another  of  this  work.  To  understand  how  to 
improve  upon  this  theory,  we  must  look  at  the  physics  of  the  problem. 

We  begin  by  taking  a  close  look  at  the  deformation  field  of  composite  plates.  It  is 
assumed  that  the  layers  of  the  laminate  are  perfectly  bonded  together  so  that  no  slip  can 
occur  between  them.  The  deformation  of  a  normal  to  the  midplane  of  the  laminate  is 
no  longer  a  smooth  continuous  function  as  it  was  for  an  isotropic  plate  as  depicted  in 
Figure  1.2,  but  instead  is  piecewise  continuous  (not  piecewise  linear,  however).  This  can 
easily  be  understood  by  considering  that  the  material  properties  can  change  drastically  at 
the  layer  interfaces.  One  can  gain  better  understanding  into  the  physics  of  the  problem  by 
studying  it  in  such  manner.  Next,  we  let  iq  be  the  three  displacements  at  any  point,  each 
being  a  general  function  of  xj.  (Here  i,j  ~  1,2,3  and  Xj  represents  the  three  coordinate 
directions.)  The  deformation  field  depicted  in  Figure  1.4  can  be  deduced  by  the  following 
items: 

1.  The  functions  Ui(xj)  must  be  continuous  in  xj  to  satisfy  compatibility. 

2.  In  the  1-2  plane,  since  there  is  no  slip  or  gaps  between  layers,  the  material  on  each 
side  of  the  interface  must  have  the  same  displacement,  and  hence,  displacement 
gradient.  In  other  words,  u,iQ  (here  a  —  1,2)  must  be  continuous  across  the  layer 
interface.  Let  us  write  this  as  ufa  =  u~Q. 
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3.  Since  ea0  =  \(ua%0  +  uf ja),  then  2  above  implies  ea0  must  be  continuous  across  the 
interface,  ie.  c+j  =  e~0. 

4.  Since  =  ea(3'  °ap  £  ca0  *n  genera^-  This  is  due  to  the  possibility  of  the 
material  constants  changing  across  the  layer  interface. 

5.  Based  simply  on  Newton’s  Third  Law,  a$j  must  be  continuous  across  the  layer 
interface,  ie.  =  a3j- 

6.  From  item  5  we  find  e3-  ^  c3j.  This  again  is  due  to  the  changing  of  material 
properties  from  layer  to  layer. 

7.  Since  ej,  ^  e3-,  and  €3 j  =  ^(u^j  +  Uj, 3),  then  u£3  ^  u~3.  (See  also  item  2  above.) 

These  seven  ideas  are  very  important  when  studying  composite  laminates.  If  one  uses  a 
displacement  based  approach  in  analyzing  laminated  plates  these  ideas  should  be  kept  in 
mind  as  the  model  is  developed.  In  doing  so,  the  limitations  and  strengths  of  any  specific 
theory  can  be  realized.  We  can  summarize  the  above  items  into  one  general  observation 
worthy  of  special  note.  This  is  as  follows: 

Observation  1  The  displacement  field,  Ui(xj),  for  a  laminated  composite  plate  must  be 
continuous  for  all  xj,  but  it  need  not  have  a  unique  u,  z  across  the  interfaces  of  the  lamina. 

This  observation  tells  us  in  no  uncertain  terms  what  form  our  displacement  field  takes.  An 
exact  solution  to  the  laminated  composite  plate  cannot  be  achieved  if  this  observation  is 
violated.  This  observation  is  supported  by  exact,  three-dimensional  analysis  of  composite 
laminates  published  by  Pagano  (1969)  [96]  and  (1970)  [97]. 
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1.5  Evaluation  of  Current  Analysis  Methodologies 

1.5.1  Higher  order  displacement  fields 

As  discussed  in  Section  1.3  there  have  been  many  attempts  at  improving  the  analysis  of 
laminated  composite  plates.  A  large  number  of  these  concentrated  upon  trying  to  gain 
more  accuracy  by  carrying  more  terms  in  the  series  expansion  of  the  displacements.  In 
other  words,  higher  order  SDPT’s  are  used  as  in  eqns  (  1.15).  One  of  the  disadvantages  of 
the  higher  order  approach  is  that  the  number  of  unknowns  in  the  problem  quickly  becomes 
large.  In  addition,  it  becomes  difficult  to  physically  understand  and  to  prescribe  boundary 
conditions  for  these  additional  terms.  The  order  of  the  expansion  by  no  means  needs  to 
be  limited  to  three  as  is  shown  in  eqns  (  1.15).  A  few  of  the  works  in  this  area  have  been 
published  by  Nelson  and  Lorch  (1974)  [84],  Kant  et  al  (1988)  [46],  Murty  and  Vellaichamy 
(1987)  [83],  Pandya  and  Kant  (1988)  [104],  Lo  et  al  (1977)  [67,  68],  Mallikarjuna  and  Kant 
(1989)  [71,  45],  Mottram  (1989)  [78]  and  Murty  (1977)  [80], 

Modifications  of  the  higher  order  approach  are  also  quite  common  in  current  research. 
They  are  referred  to  as  the  simplified  higher  order  methods.  Through  setting  the  transverse 
stresses  on  the  top  and  bottom  surfaces  of  the  plate  equal  to  zero,  the  number  of  unknowns 
in  the  displacement  field  can  be  reduced  by  four  terms.  For  example,  the  most  common 
simplified  higher  order  method  found  in  the  literature  begins  with  cubic  forms  for  u  and  v 
and  a  constant  w  (with  respect  to  z).  This  nine  parameter  displacement  field  can  then  be 
reduced  to  one  of  only  five  parameters,  the  same  number  as  found  in  the  first  order  theory. 
This  specific  example  will  be  discussed  in  more  detail  in  Section  3.1.  Some  examples  of 
simplified  higher  order  approaches  can  be  found  in  references  published  by  Reddy  (1984) 
[113,  112],  Kant  and  Pandya  (1988)  [44],  Senthilnathan  et  al  (1988)  [124,  65],  Reddy  and 
Phan  (1985)  [117],  Murty  (1987)  [81]  and  Khdeir  (1989)  [53]. 

These  higher  order  approaches  theoretically  make  sense,  but  outside  of  the  realm  of 
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homogeneous  materials,  carrying  a  finite  number  of  higher  order  terms  does  not  coincide 
with  the  physics  of  the  problem.  A  higher  order  SDPT  has  unique  values  of  u, 3  across 
the  lamina  interfaces.  Thus,  a  higher  order  SDPT  does  not  have  the  freedom  to  fully 
comply  with  Observation  1.  The  advantage  to  a  higher  order  theory  is  that  a  traction 
free  condition  on  the  top  and  bottom  surfaces  can  be  satisfied,  and  the  need  for  a  shear 
correction  coefficient  is  reduced.  However,  accurate  prediction  of  the  transverse  stresses 
cannot  be  done  directly.  We  use  the  word  directly  here  because  the  transverse  stresses  can 
be  found  by  integrating  the  equations  of  equilibrium  through  the  thickness  of  the  plate. 
This  is  done  utilizing  the  calculated  in-plane  stresses  which  are  generally  accurate.  This 
technique  is  used  quite  often  when  transverse  stresses  are  desired.  Some  examples  of  this 
are  included  in  the  works  by  Nishioka  and  Atluri  (1979)  [85],  Murty  (1987)  [82],  Kant  and 
Pandya  (1988)  [44],  Noor  and  Burton  (1989)  [88]  and  Reddy  et  al  (1989)  [116].  Also,  note 
the  above  emphasis  on  the  word  reduced  with  regards  to  the  shear  correction  coefficient. 
Since  the  assumed  deformation  field  is  not  exact,  then  the  calculated  transverse  strain 
energy  is  going  to  deviate  from  the  true  energy,  and  a  shear  correction  coefficient  would 
be  beneficial  to  adjust  this  difference.  If  the  assumed  deformation  closely  approximates 
the  true  one,  then  the  shear  correction  coefficient  would  be  insignificant.  In  other  words, 
in  studying  an  isotropic  plate  one  would  expect  a  third  order  theory  to  give  excellent 
results  without  a  shear  correction  factor,  as  such  a  field  follows  an  elasticity  solution. 
However,  depending  upon  the  makeup  of  a  composite  laminate,  a  third  order  theory 
cannot  accurately  describe  a  piecewise  continuous  displacement  field.  Based  upon  this 
discussion,  it  is  reasonable  to  assume  that  the  first  order  SDPT  with  the  correct  shear 
correction  coefficients  can  give  better  results  for  certain  cases  than  a  higher  order  theory 
without  any  correction.  This  time,  notice  that  the  emphasis  is  placed  upon  the  word  correct 
with  regard  to  the  shear  correction  coefficients.  One  of  the  purposes  of  this  work  is  to 
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demonstrate  that  the  shear  correction  coefficients  can  be  calculated  accurately,  allowing 
the  first  order  SDPT  to  provide  as  good,  if  not  better,  results  than  the  higher  order 
approaches  with  their  added  complexities. 

1.5.2  Discrete  layer  approach 

Up  to  this  point  we  have  not  considered  the  obvious  methodology  for  the  laminated  com¬ 
posite  plate  problem.  This,  of  course,  is  to  assume  a  piecewise  continuous  displacement 
field  through  the  thickness  of  the  laminate  as  depicted  in  Figure  1.4.  Within  each  layer 
the  displacements  can  be  chosen  to  be  linear  or  of  a  higher  order  form.  This  form  of 
displacement  field  should  provide  us  with  more  accurate  results  than  a  finite  term  higher 
order  approach  and  is  often  referred  to  as  a  discrete  layer  approach.  The  discrete  layer 
approach  is  nothing  more  than  modeling  each  individual  layer  of  the  laminate  as  a  sepa¬ 
rate  plate.  The  obvious  drawback  to  such  an  approach  is  that  the  number  of  unknowns  in 
the  problem  becomes  tied  to  the  number  of  layers  in  the  laminate.  A  problem  can  quickly 
become  intractable  for  thick  plates  with  a  large  number  of  layers.  A  few  of  the  works 
published  in  this  area  include  those  by  Srinivas  (1973)  [130],  Reddy  et  al  (1989)  [116], 
Barbero  and  Reddy  (1990)  [6]  and  Alam  and  Asnani  (1984)  [2,  3]. 

Just  as  there  was  a  simplified  version  of  the  higher  order  approach,  there  is  also  a 
simplified  discrete  layer  approach.  In  assuming  a  piecewise  continuous  displacement  field, 
the  number  of  unknowns  can  be  reduced  by  enforcing  displacement  continuity,  as  well  as 
transverse  shear  traction  continuity,  at  the  layer  interfaces.  In  this  manner  the  number 
of  unknowns  can  be  made  to  be  independent  of  the  number  of  layers.  (This  will  be 
developed  in  more  detail  in  Section  3.2.2.)  The  simplified  discrete  layer  approach  has 
been  demonstrated  to  provide  accurate  results  for  thick  composite  laminates  and  has  the 
potential  to  become  an  efficient  and  powerful  tool  in  the  analysis  of  laminated  composite 
plates.  A  few  examples  of  the  simplified  discrete  layer  approach  include  Sciuva  (1987) 
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[123],  Mawenya  and  Davies  (1974)  [73]  and  Lee  et  al  (1990)  [59]. 

It  is  interesting  to  note  that  a  piecewise  continuous  displacement  field  results  in  a 
smooth  transverse  stress  field  rather  than  a  piecewise  continuous  one  through  the  thick¬ 
ness.  This  is  because  when  the  interlaminar  transverse  stress  continuity  is  enforced  across 
a  lamina  interface,  either  in  the  discrete  layer  or  simplified  discrete  layer  theory,  the 
derivatives  (slopes)  with  respect  to  z  of  these  functions  are  also  the  same.  Thus,  a  smooth 
function  through  the  thickness  results.  The  end  result  is  that  in  order  to  calculate  the 
transverse  stresses,  one  must  rely  upon  integrating  the  equations  of  equilibrium  as  dis¬ 
cussed  in  Section  1.5.1. 

1.5.3  Finite  Element  Analysis  of  Composite  Materials 

Finite  element  analysis  of  laminated  plates  began  with  the  work  of  Pryor  and  Barker 
(1971)  [42],  In  their  work  they  developed  an  element  based  upon  the  deformation  field 
given  in  eqns  (  1.2)-(  1.4),  and  established  seven  degrees  of  freedom  at  each  node  (ua,  va, 
wa,  4>x,  <f> y,  7 xz  and  7 y2).  This  basic  development  has  been  the  fundamental  basis  for  much 
of  the  finite  element  work  done  with  composite  materials  to  date,  with  the  only  difference 
being  that  five  degrees  of  freedom  per  node  are  generally  used  (ua,  v0,  wa,  <f>z  and  <f>y). 
This  basic  approach  is  the  back-bone  of  composite  analyses,  and  most  commercial  finite 
element  codes  are  based  upon  this  theory.  Many  different  elements  and  implementations 
have  been  devised,  but  generally  the  first  order  shear  deformation  theory  is  the  basis 
behind  them.  Some  of  the  works  published  in  the  past  several  years  using  a  first  order 
shear  deformation  approach  include:  Moser  and  Schmid  (1989)  [77],  Craig  and  Dawe 
(1987)  [24],  Kumar  and  Rao  (1988)  [57],  Fuehne  and  Engblom  (1989)  [26],  Lardeur  and 
Batoz  (1989)  [58],  Zienkiewicz  and  Lefebvre  (1988)  [147],  Onate  and  Suarez  (1983)  [95], 
Reddy  (1980)  [110],  Suresh  et  al  (1979)  [103]  and  Irons  and  Zienkiewicz  (1970)  [38].  As 
one  might  expect,  the  first  order  theory  is  burdened  with  the  shear  correction  coefficient 
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problem  as  discussed  in  Section  1.3.  This  problem  is  even  more  difficult  with  laminates 
than  with  homogeneous  plates,  because  the  coefficients  not  only  depend  upon  geometry, 
but  also  on  ply  orientation  and  material  properties.  In  addition,  because  of  this,  different 
coefficients  are  also  needed  for  different  directions  in  the  plate.  However,  despite  this 
drawback,  the  first  order  approach  has  remained  popular  because  of  its  low  computational 
cost.  Implementation  of  higher  order  approaches  has  also  been  well  documented  in  the 
published  work.  These  include  works  by  Kant  et  al  (1988)  [46,  44,  104],  Mallikarjuna  and 
Kant  (1989)  [71]  and  Reddy  (1989)  [115].  The  advantage  of  these,  as  discussed  previously, 
is  the  reduced  need  for  shear  correction  coefficients,  while  the  disadvantage  is  the  increased 
complexity.  The  simplified  higher  order  approaches,  even  with  lower  numbers  of  degrees  of 
freedom,  have  increased  complexity  in  that  they  have  increased  continuity  requirements. 
This  will  be  discussed  in  more  detail  in  Section  3.2. 

The  discrete  layer  and  simplified  discrete  layer  approaches  have  also  been  employed 
in  finite  element  analyses.  The  full  discrete  layer  approaches  are  computationally  more 
intense  than  other  methods  and  are  not  widely  found  in  the  literature.  Examples  include 
those  by  Reddy  et  al  (1989)  [116]  and  Barbero  and  Reddy  (1990)  [6].  Examples  of  the 
simplified  discrete  layer  approaches  are  computationally  efficient  and  include:  Lee  et  al 
(1990)  [59]. 

All  of  the  above  mentioned  finite  element  analyses  are  displacement  based  approaches. 
Along  with  these  come  the  lack  of  direct  transverse  stress  calculations  and  increased 
continuity  requirements  for  the  finite  element  model.  These  inherent  problems  have  led  to 
much  research  in  the  areas  of  mixed  and  hybrid  methods  which  eliminate  these  problems. 
These  methods  have  demonstrated  excellent  results,  but  again  at  the  cost  of  excessive 
computations.  Research  in  these  areas  include:  Putcha  and  Reddy  (1986)  [106],  Jing  and 
Liao  (1989)  [39],  Spilker  (1982)  [128],  Spilker  et  al  (1977)  [129],  Mau,  Tong  and  Pian 
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(1972)  [72]  and  Liou  and  Sun  (1987)  [66], 

As  a  final  note,  there  have  been  a  few  survey  and  overview  articles  by  Reddy  (1985) 
[114],  (1981)  [111],  (1989)  [115],  which  have  appeared  in  the  past  several  years.  The 
interested  reader  may  find  them  beneficial  if  more  information  is  desired. 

1.6  Analysis  Approach 

1.6.1  Improved  First  Order  SDPT 

As  discussed  above  in  Section  1.5.1  the  first  order  SDPT  has  the  potential  to  provide 
accurate  results  as  long  as  the  correct  shear  correction  coefficients  earn  be  calculated.  In 
Chapter  II  we  will  develop  a  finite  element  method  based  upon  first  order  SDPT  which  can 
accurately  calculate  the  vibration  characteristics  of  thick  laminated  composite  plates.  The 
method  will  include  the  calculation  of  accurate  shear  correction  coefficients  by  comparing 
the  first  order  transverse  strain  energy  to  a  more  accurate  strain  energy.  This  more 
accurate  strain  energy  will  be  calculated  based  upon  the  transverse  stresses  found  from 
integrating  the  equations  of  equilibrium  through  the  thickness  of  the  plate. 

1.6.2  Simplified  Discrete  Layer  Approach  with  a  Least  Squares  Element 

Chapter  III  of  this  work  will  develop  two  new  Least  Squares  finite  elements.  The  elements 
will  utilize  a  unique  Least  Squares  method  to  allow  an  element’s  displacement  functions 
to  approximate  Cl  continuous  functions  on  the  boundaries  of  the  element.  This  technique 
will  allow  the  element  to  behave  as  one  with  C1  continuity.  Thus,  it  will  enable  the  use  of 
displacement  functions  which  normally  would  require  C 1  continuous  interpolation  func¬ 
tions.  The  Least  Squares  element  would  be  applicable  to  either  a  simplified  higher  order 
approach  or  a  simplified  discrete  layer  approach  with  a  simplified  higher  order  piecewise 
continuous  function  as  the  basis  for  the  displacement  field.  These  forms  of  displacement 
functions  will  be  shown  to  contain  differentials  of  the  out-of-plane  displacement,  which  up 
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until  now  could  not  easily  be  analyzed  without  the  use  of  C1  continuous  elements.  The 
simplified  discrete  layer  field,  based  upon  the  simplified  higher  order  approach,  will  be 
used  to  demonstrate  the  technique. 

1.6.3  Effect  of  Fiber  Orientation  and  Stacking  Sequence  on  Fundamental 
Frequency  and  Stability 

The  Least  Squares  finite  element  technique  will  be  used  to  perform  a  study  of  the  effects 
of  fiber  orientation  angle,  stacking  sequence,  boundary  conditions,  aspect  ratios  and  pre¬ 
stressing  on  the  fundamental  frequencies  and  buckling  loads  of  composite  laminates.  The 
study  will  concentrate  upon  optimizing  frequency  and  buckling  load  while  varying  the 
other  parameters.  The  purpose  of  the  study  is  to  provide  new  information  for  the  design 
engineer  to  aid  in  the  designing  of  composite  laminated  materials. 


CHAPTER  II 


AN  IMPROVED  FIRST  ORDER  SHEAR 
DEFORMATION  THEORY  THROUGH  THE 
USE  OF  THE  PREDICTOR  CORRECTOR 

TECHNIQUE 

2.1  Analysis  Overview 

In  this  part  of  the  work  we  will  investigate  the  feasibility  of  developing  a  finite  element 
model  which  is  based  upon  a  first  order  SDPT.  The  method  stands  out  from  any  other  finite 
element  models  of  this  type  because  it  will  include  the  ability  to  calculate  accurate  shear 
correction  coefficients  for  any  particular  laminate  being  considered.  With  accurate  shear 
correction  coefficients,  the  first  order  SDPT  will  give  very  good  results.  This  technique 
has  been  successfully  applied  analytically  by  Noor  and  Peters  (1989)  [91]  and  Noor  and 
Burton  (1989)  [88],  (1990)  [89],  but  implementation  of  the  technique  into  a  finite  element 
analysis  has  not  been  published. 

In  the  following  sections  we  will  first  develop  a  first  order  plate  bending  model  for  a 
general  composite  laminate.  This  development  will  be  the  basis  for  the  development  of 
the  finite  element  model  to  be  used  for  a  vibration  analysis.  The  finite  element  model  will 
be  developed  to  include  the  integration  of  the  equilibrium  equations  through  the  thickness 
of  the  plate  to  obtain  the  transverse  stresses.  These  stresses  will  then  be  used  to  calculate 
the  transverse  strain  energy  which  will  be  compared  to  those  obtained  using  the  first  order 
theory.  This  comparison  will  yield  the  above  mentioned  shear  correction  coefficient  which 
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then  will  be  used  to  calculate  an  updated  natural  frequency  and  mode  shape. 

2.2  Preliminary  Development 

2.2.1  Fundamental  Equations 


As  presented  and  discussed  above  in  Section  1.2,  the  assumed  displacements  for  a  first 
order  SDPT  are  represented  by  eqns  (  1.2)-(  1.4).  These  equations  will  form  the  basis  for 
our  analysis.  In  terms  of  these  variables,  eqns  (  1.7)-(  1.12)  represent  the  strains. 

Next,  we  can  begin  to  introduce  some  of  the  standard  equations  use  ’  .  hen  working 
with  fiber  reinforced  composite  materials  (see  texts  by  Jones  [41]  and  Christensen  [22]. 
First  we  define: 


for  a  single  layer  of  a  composite  material  related  through  the  constitutive  equation 


{a}  =  [C]{e}  (2.1) 

where  [C]  are  the  stiffnesses  in  the  1-2-3  coordinate  system  as  depicted  in  Figure  1.3.  If 
we  assume  two  orthogonal  planes  of  material  property  symmetry  for  the  material  under 
consideration,  then  [C]  takes  the  standard  form: 


Cn  C 12  C i3  0  0  0 

C\2  C22  C23  0  0  0 

C13  C23  C33  0  0  0 

0  0  0  Ce  s  0  0 

0  0  0  0  C44  0 

0  0  0  0  0  C55 


The  values  for  the  Cij  are: 


(2.2) 


Cn  =  #1  (1  —  ”32*23)  /A,  C12  =  £2  (1^12  +  ”13*32)  (A,  C44  =  G\z 

C12  =  £3  (l/!3  +  Z/121/23)  M,  C22  =  £2(1-  ^31^13)  M,  C55  =  G23 

C23  =  Ez{v2Z  +  ”\3”2\)  j A,  C33  =  £3  (1  -  1^12^21)  /A,  C$6  =  G12 
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(2.3) 

where 

A  =  (1  —  U23U32  ~  ”12”2l  —  ”13”31  —  V\2”23”3l  —  1^13^32*^21) 

We  also  define  the  Poisson’s  ratio,  to  be  the  ratio  of  the  deformation  in  the  j  direction 
to  that  in  the  i  direction  when  pulled  in  the  i  direction. 

The  next  step  is  to  then  transform  [C]  into  the  x,  y,  z  coordinate  system  to  get  the 
form: 


axx  Cn  C\2  C13  C 16  0  0 

ayy  C\2  C22  <?23  <?26  0  0 

,  azz  (  _  <?13  <?23  C33  C36  0  0 

axy  Ci6  C26  <?36  <?66  0  0 

< Txz  0  0  0  0  C44  C45 

.  ayz  .  .0  0  0  0  C45  C55 


where  (n  =  sin  9,  m  —  cos  9): 


Cu  =  m4Cn  +  2m2n2(C12  +  2C66)  +  n4C22 

C12  =  m2n2  (Cn  +  C22  -  4Cqq)  +  (m4  4-  n4)  C12 

C22  =  n4Cn  +  2m2n2  (C12  +  2Cee)  +  m4C22 

Ci6  =  -mn  [n2Cn  -  rr^Cn  -  (m2  -  n2)  (C12  +  2C66)] 

C26  =  -mn  [m2Cn  -  n2C22  +  (m2  -  n2)  (C12  +  2C66)j 
^66  =  m2n2  (Cn  +  C22  -  2C12)  +  (m4  -  n4j  C66 


(2-5) 


24 


Ci3  =  m2Ci3  +  n2C23  C44  =  m2C44  +  n2C55 
C23  =  n*C\z  4-  Tn2C23  C36  =  (C32  —  C31)  nm 
C45  =  (C44  —  C55)  nm  C55  =  m2Czb  +  "r?C\\ 

C33  =  C33 

We  can  next  eliminate  the  azz  equation  from  eqn  (  2.4).  One  can  solve  for  ezz  from  eqn 
(  2.4)  to  get 


^zz  —  q  {&zz  Ciz€Xx  Czztyy  C36T11/) 


(2.6) 


which  can  then  be  substituted  back  into  eqn  (  2.4)  to  yield: 


Q11 

Ql2 

Ql6 

0 

0  ' 

’  ’ 

crr 

'  ?13  ' 

ayy 

Q 12 

Q22 

Q26 

0 

0 

eyy 

ST 

C23 

17  zy 

&xZ 

>  — 

Ql6 

0 

Q26 

0 

C?66 

0 

0 

Q44 

0 

Q45 

< 

Try 

Txz 

U  zz 

C33 

C36 

0 

.  ay*  , 

0 

0 

0 

Q45 

Qss  . 

.  Tyz  . 

.  0  . 

(2.7) 


The  values,  Q,j,  are  the  reduced  stiffness  coefficients  developed  by  Whitney  (1969)  [141], 
and  are  defined  in  terms  of  the  C,j  as: 


Qij  =  Cij  -  ( Ciz  ■  Cjz)  /C33  i,j  =  1, 2, 6 
Qij  —  Cij  i,j  =  4,5 


(2.8) 


Next,  before  moving  on,  it  will  be  convenient  to  have  [ Q ]  partitioned  into  two  matrices. 


[Q]  and  [Q],  defined  as 


[Q] 


[Q]  0 

0  [Q] 


(2.9) 
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The  reason  for  breaking  [Q]  up  like  this  will  be  apparent  in  the  next  section.  Lastly,  from 
this  point  on,  we  will  drop  the  term  containing  azz  in  eqn  (  2.7).  In  other  words  we  assume 
that  the  transverse  normal  stress  has  little  effect  on  our  solution. 

2.2.2  Internal  Forces  and  Moments 


At  this  point  it  is  convenient  to  define  the  internal  force  and  moment  relations  for  a 
composite  laminate.  The  forces  per  unit  width  of  the  plate  are  found  by  integrating  the 
stresses  through  the  thickness  to  get 


where  we  have  introduced  k{  and  kj,  the  shear  correction  coefficients,  which  will  be  dis¬ 
cussed  in  more  detail  later.  Similarly,  the  moments  per  unit  width  are 


\  1 

fhl  2 

'  <7„  ' 

rh/2 

=  m  =  /  J 

_ 

°yy 

o: 

ii 

N 

■« 

[  MZy  t 

J-h/2 

j 

J-h/2 

.  . 

Next,  we  use  eqns  (  1.7)-(  1.12)  to  express  the  strains  in  terms  of  midplane  strains  and 
curvatures.  The  result  is 


e°  +  zKx  ) 
ty  +  ZKy  f 

To  "h  J 


dz=[  ([<?]  {e0}  +  [Q\z{K})dz 

J-h/2 


(2.13) 
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(Q) = (**)*  f_'jQ\ { X% } * = <**>* £'>1 w *  (2.i4) 

rh/2  t°x  +  ZKx  'j  h/2 

{M}  -  J  h  \Q\  '  €y  +  ZKy  \  z  dz  =  J  h  2  ([<?]Z  +  ^Z  (2-15) 

7o  +  2Kiy  J 

In  the  above  three  equations,  the  [Q]  and  [Q]  are  piecewise  constant  with  respect  to  the 
integration  through  the  thickness  of  the  laminate.  Thus,  as  is  commonly  done,  we  establish 
matrices  [j4],  [j4],  [ B\  and  [£>]  whose  elements  are  defined  by: 


rh/2  » 

■A-ij  =  /  Qijdz  —  J^(Qij)k  '  4 

J-h/2 


Aij  =  (*,*,)’  jh/2  Qijdz  =  {kik])hj^{QlJ)k-tk 

J~h /2  it=l 

/•/j/2  » 

Aj  =  /  Qijz  '  =  y~l(Q»7^fc  ‘  tk  '  Zk 

J-h/2  fc=i 

/■h/2  2  /  f3  \ 

Aj  =  2  Q'iz  '  =  +  ^2  J 


(2.16) 

(2.17) 

(2.18) 

(2.19) 


where  zk  =  zk  -  zk-\,  and  4  and  z*.  are  defined  as  in  Figure  2.1.  Finally,  since  {f0} ,  {«} 
and  {t?}  are  independent  of  2,  they  can  be  pulled  outside  of  the  integral.  The  end  result 
is  a  convenient  representation  of  the  internal  forces  and  moments  in  the  laminate: 

N  |  r  A  B  0  ]  f  fo  ' 

M  -  =  D  D  0  i  k  •  (2.20) 

Q  0  0  A  I  0 

V  L  J  \ 
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Z 

A 


Layer  n 
Layer  n -1 


Figure  2.1:  Laminate  Geometry 


These  equations  will  now  be  used  in  the  finite  element  formulation. 

2.3  Finite  Element  Formulation 

2.3.1  Element  Mathematical  Considerations 

As  discussed  in  the  previous  section,  the  basis  for  this  analysis  is  a  first  order  shear 
deformation  theory.  The  displacements  through  the  thickness  of  the  plate  will  be  described 
with  three  displacements  and  two  rotations.  As  a  result,  a  two  dimensional  element  is  all 
that  is  required  rather  than  a  three  dimensional  solid  one.  The  quadratic,  eight  noded, 
isoparametric  element  is  often  used  in  modeling  plates  with  accurate  results.  For  this 
reason  it  will  be  used  for  this  analysis.  The  eight  noded  element  will  have  the  three 
displacements  and  two  rotations  at  each  node  for  a  total  of  40  degrees  of  freedom.  The 
element  is  shown  in  Figure  2.2. 

The  element  parameters  will  be  described  through  standard  shape  functions  derived 
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* 


Figure  2.2:  Eight  Noded  Isoparametric  Element 

for  the  serendipity  family  of  elements  as  given  in  any  standard  finite  element  text.  (See 
Zienkiewicz  [146],  Cook  [23],  Tong  and  Rossettos  [136]  or  Bathe  [9].)  The  eight  shape 
functions  are 

Ni  =  \  (1  +  Wi)  (1  +  «i)  («.  +  VVi  ~  1)  *  =  1, 2, 3, 4 
Ni  =  ±(  l  +  w)  (l  -  Z2)  *  =  5,7  (2.21) 

=  U  l-i?2)(l  +  ^I)  *  =  6,8 

where  here  rji  and  Eire  the  coordinates  of  the  ith  node,  Eind  T)  and  £  are  the  coordinates 
of  any  desired  point  in  the  element.  The  diseidvantage  of  this  element,  which  will  be 
discussed  in  detail  later,  is  that  it  is  only  capable  of  having  a  quadratic  variation  of 
the  primEuy  variables  throughout  it.  Thus,  the  stresses  will  be  linear  and  derivatives 
of  stresses  constant.  This  restriction  can  be  overcome,  but  as  mentioned  above,  will  be 


discussed  later. 
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Finding  an  extremum  of  this  functional  leads  to  the  stiffness  and  mass  matrices. 

2.3.2  Strain  energy  development 

In  the  above  discussion,  eqn  (  2.22)  can  readily  be  adapted  for  the  case  of  the  composite 
laminate.  If  we  perform  an  integration  through  the  thickness  of  the  laminate,  just  as 
performed  in  the  development  of  eqn  (  2.20),  it  can  be  shown  that  the  expression  for  the 
strain  energy  of  the  system  becomes 

K— £/.({£  H  « }  +  {<wTw)«  <2-25> 

where  it  is  now  written  in  terms  of  the  internal  forces  and  moments  and  is  an  area  integral 
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rather  than  a  volume  integral.  The  goal  is  to  write  this  expression  ultimately  in  terms  of 
the  nodal  degrees  of  freedom  of  the  element.  Prom  eqn  (  2.20)  {N}  ,  {M},  and  {Q}  can 
be  written  as 


'  A  B 
B  D 


{Q}  = 


(2.26) 

(2.27) 


Next,  recalling  the  definitions  of  {€°},{k},  and  {t?},  (see  eqns  (  1.7)-(  1.12)  ),  we  can 
express  them  as 


uo,x 

'  1 

0 

0 

0 

0 

0 

0 

0  ' 

vo,y 

0 

0 

0 

1 

0 

0 

0 

0 

vo,x  +  uo,y 

0 

1 

1 

0 

0 

0 

0 

0 

i 

$1,1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

.  $x,y  +  <j>yyx  , 

.  0 

0 

0 

0 

0 

1 

1 

0  . 

uo,x 
uo,y 
vo,  x 
vo,y 
4>x,x 
<f>x,y 
<f>y,x 


(2.28) 


{*) 


W,x  +4>x 

w,y  ~^~4>y 


} 


10  10 
0  10  1 


4>x 

4>y 


w 


o,y 


(2.29) 


Next,  and  j,  which  are  expressed  in  terms  of  the  derivatives  with  respect  to  x 

and  y,  need  to  be  written  in  terms  of  the  derivatives  with  respect  to  7;  and  £.  This  can 
easily  be  done  through  the  standard  use  of  the  Jacobian  matrix,  [J],  which  relates  the 
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derivatives  of  the  two  coordinate  systems.  If  we  define 

m  =  [j]-1 

then  we  can  write: 

r  -ft*.,) 

r  o  o  o 

fc}  =  o  o  r  o  I  ?  l  =  [//)Kd  (2-3°) 

o  o  o  r  J  : 

4>x 

<t>y 

w,n 

We 

where  I  used  in  |//j  is  the  identity  matrix.  The  final  step  is  to  now  represent  j  and 

1 6'^  |  in  terms  of  the  forty  nodal  degrees  of  freedom.  First,  we  establish  the  order  of  the 

nodal  degrees  of  freedom  (eight  sets  of  the  five  displacements)  to  be 

(«o)l  ' 

(*»<>)  1 
(»o)l 

{A}  =  (Ml  (2.32) 

(<f>y)  1 
.  (<f>y) 8  . 


so  that  through  the  use  of  the  shape  functions  we  can  represent 


'  Nhv  0  0  0  0  N2,r,  0  •••  0 

Nu  0  0  0  0  N2<z  0  •••  0 

o’  Nhr]  0  0  0  o’  N2,r,  0 

0  iVu  0  0  0  0  N2£  •••  0 

M  -  ;  ; 


0  0  0  0  Ni4  •••  •••  Ns 
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[777]  {A} 


'  0 

0 

0 

Ni 

0 

0 

0 

0 

n2  • 

•  0  ‘ 

0 

0 

0 

0 

Ni 

0 

0 

0 

0  • 

•  0 

i 

1  A  1 

0 

0 

0 

0 

0 

0 

*2,, 

0  • 

•  0 

1 

0 

0 

*1* 

0 

0 

0 

0 

*24 

0  • 

0 

l  J 

(2.33) 


=  [//i]{  A} 


(2.34) 


So  now  we  have  in  toted 


Co 

K 


=  [7] [77] [777]  {A}  =  \J3\  {A} 


(2.35) 


{1?}  =  [7][7-7][777]{A}  =  [/0]{A} 


(2.36) 


Substituting  these  into  eqns  (  2.26)-(  2.28),  the  equations  for  the  internal  forces  and 
moments  we  get: 


r  n  i 

'  A 

B 

lM  J 

B 

D 

Ifl  {A} 


(2.37) 


{Q }  =  [a]  &9]{A} 


(2.38) 


With  these  representations,  the  expression  for  total  strain  energy  can  now  be  expressed  as 


U, 


-i  /.( 


{a  yip? 


A  B 
B  D 


1  T 


[P\  {A}  +  {A}7  ffi[A)T\j3)  {A}  I  dA  (2.39) 
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2.3.3  Kinetic  energy  development 

The  mass  matrix  for  the  dynamic  analysis  of  the  composite  plate  comes  from  the  expression 
for  the  kinetic  energy  of  the  plate  given  in  eqn  (  2.23)  .  The  time  derivatives  of  the 
displacements  in  eqns  (  1.2)-(  1.4)  become 


u  =  u0  +  z4>x 

v  =  v0  +  z4>y  (2.40) 


w  =  wa 


Substituting  these  expressions  into  eqn  (  2.23)  gives 


Uke~\JvP  +  "o  +  2z  +  Votv)  +^l  +  z2  (4>2X  +  dV  (2.41) 


If  the  displacement  functions  cure  assumed  to  be  harmonic  functions  of  time,  each  one  can 
be  assumed  to  be  premultiplied  by  e'ut.  With  this,  the  time  derivatives  can  be  written 


as  the  displacements  themselves  multiplied  by  iu.  To  aid  in  the  development  of  the  mass 


matrix  we  define  the  following  expressions: 
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.  0 
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1 

{6}  =  iu,[Il}{6} 


(2.42) 


{6}=iu;[I2}{6} 


(2.43) 
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'00010' 

0  0  0  0  1 

iu  0  0  0  0  0  {6}  =  iu  [/3]  {6}  (2.44) 

1  0  0  0  0 

.0  1  0  0  0. 

'1  0  0  0  0' 

0  10  0  0 

iu  0  0  0  0  0  {6}  =  iu  [/4]  {6}  (2.45) 

0  0  0  1  0 

.0  0  0  0  1. 


Now,  with  these  equations  defined,  we  can  write  eqn  (  2.41)  as 


Uk,  =  jvp({i)T  [h]T  [/,]{«} 

+z2  {«}T  (//'  [/2]  {«}  +  z  {«(r  [/3f  [/„]  («})  dV  (2.46) 

which  simplifies  to,  (due  to  the  nature  of  the  matrices), 


Uke  =  -\“2JvP  ({*}T  [A]  {6}  +  22  {^}r  [h]  {6}  +  Z  {6}t  lh)  {6})  dV  (2.47) 


In  this  equation,  only  p  and  z  vary  through  the  thickness,  so  as  before,  we  can  easily 
perform  the  integration  in  the  z  direction,  and  reduce  it  down  to  an  area  integral.  The 
result  is 


Uke  =  -^2  JA  (p  {Sf  [/i]  {6}  +  I  {£}T  [/2]  {6}  +  R  {6}t  lh\  {«})  dV  (2.48) 

where 

rh/2  n 

P  =  /  pdz  =  YJPk-tk  (2.49) 

J~W  1 1 
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R 


I 


pkz  -dz=Y^Pk-tk-Zk 
k=l 


pz2  ■ dz  =  Ylpk 

k= 1 


(“ 


4 + 


(2.50) 


(2.51) 


The  final  form  for  Uke  follows  after  writing  {5}  in  terms  of  {A}  and  the  shape  functions. 
We  have 
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{A} 


=  [N){A} 


Now  we  can  write  our  final  form  of  Uke  as 


(2.52) 


Uke 


P  0 
0  I 


[iV]{A} 


+  {A}T[iVlT 


0  R 
R  0 


[N]  {A}  dA 


(2.53) 


where  P  [Ji]  and  /  [fj]  have  been  combined  into  one  matrix.  With  this  equation,  along 
with  eqn  (  2.39),  we  have  what  is  needed  to  minimize  the  energy  functional  and  get  an 
expression  for  the  stiffness  and  mass  matrices. 

2.3.4  Element  stiffness  and  mass  matrices 

In  the  last  two  sections  we  developed  expressions  for  both  the  strain  energy,  eqn  (  2.39)  , 
and  kinetic  energy,  eqn  (  2.53)  ,  for  a  general  composite  laminate.  Substituting  these  two 
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•  equations  into  eqn  (  2.24)  ,  the  expression  for  total  energy,  results  in 


n  =  \ fA({*)T [fif  g  bd  i^HA}  +  (A}r[4]Wl^]{A}j^ 

-r2L{iAf{N]T[ o  ?]IW|<A> 

+  {A}r[JV]T  °  *  [JV]  f  A}  j  dA  (2.54) 

If  we  extremize  this  equation  with  respect  to  the  nodal  degrees  of  freedom,  the  result  is 

/([0]T  B  D  K  +  WW 

“"’/rfM  0  ?]<"] 

+[N)T  “  J  [JV]j{A)^=0  (2.55) 

This  equation  is  of  the  form 

([A;]  -  w2  [m])  {A}  =  0  (2.56) 

giving  us  our  eigenvalue  problem  to  solve  for  the  natural  frequencies  and  mode  shapes  of 
our  problem.  From  eqn  (  2.55)  the  elemental  stiffness  matrix  is 

[*]  =  L  ((/?]T  B  D  l P}  +  \P}TlA}T[P}^JdA  (2.57) 


and  the  elemental  mass  matrix  is 
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H  = 


P  0 
0  I 


[N]  +  [N)t 


0  R 
R  0 


[N^j  dA 


(2.58) 


At  this  point  it  is  important  to  remember  that  the  integration  must  be  done  in  the  rj  - 
£  coordinate  system  since  the  displacements  have  been  represented  in  terms  of  the  shape 
functions.  Thus,  the  differential  area  element  becomes: 

dA  =  ||J||dj7<i£ 

2.4  Analysis  Method 

2.4.1  Global  matrix  development 

In  Section  2.3.4  we  developed  expressions  for  the  elemental  stiffness  and  mass  matrices  for 
a  laminated  composite  plate  which  was  based  upon  an  eight  noded  isoparametric  element. 
The  expressions,  given  by  eqn  (  2.57)  and  eqn  (  2.58),  are  the  basis  for  developing  the 
global  system  of  equations  in  the  finite  element  problem. 

The  process  of  building  global  stiffness  and  mass  matrices  follows  the  standard  proce¬ 
dure  of  discretizing  the  domain  in  question  to  establish  an  elemental  mesh.  The  elements 
and  nodes  are  numbered,  and  the  integration  of  eqn  (  2.57)  and  eqn  (  2.58)  is  performed 
over  each  element.  Thus,  we  build  the  element’s  stiffness  and  mass  matrices.  The  results 
for  each  element  are  then  assembled  into  a  globed  mass  and  stiffness  matrix  resulting  in  a 
matrix  equation  of  the  form: 


([/?]-o,2[mi){a}  =  o. 


(2.59) 


In  carrying  out  the  integration  of  eqn  (  2.57)  and  eqn  (  2.58)  one  cem  easily  employ  the 
Gauss  queidrature  technique.  Unfortunately,  the  number  of  integration  points  used  for  the 
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different  terms  in  the  equations  varies.  For  example,  in  eqn  (  2.58)  the  integrand  contains 
the  squares  of  the  shape  functions  given  in  eqn  (  2.21).  The  result  is  an  expression  with 
terms  on  the  order  of  q4  and  £4.  Thus,  to  get  an  accurate  integration  of  the  mass  and  inertia 
terms,  a  three  by  three  point  Guassian  integration  must  be  performed.  Unfortunately,  eqn 
(  2.57)  presents  another  case.  The  first  terms,  which  represent  the  bending  stiffnesses,  are 
on  the  order  of  172  and  £2,  so  can  be  integrated  with  a  two  by  two  scheme.  The  next 
expression  contains  the  shear  stiffness  terms  which  are  on  the  order  of  r/4  and  £4,  as 
well  as  t/2  and  £2.  The  fourth  order  terms  result  from  the  fact  that  the  transverse  shear 
contains  <f>x  and  <f>y  and  not  their  derivatives.  However,  despite  the  higher  order  of  the 
shear  stiffness  terms,  it  is  often  recommended  that  a  one  point  integration  scheme  be  used 
to  give  the  best  results  when  solving  eqn  (  2.59)  for  the  natural  frequencies.  This  need  for 
under  integration  follows  from  the  presence  of  the  shear  locking  phenomenon  as  described 
in  the  text  by  Hughes  [36].  When  the  shear  terms  are  integrated  with  a  lower  order  than 
the  bending  terms,  it  is  referred  to  as  reduced  selective  integration.  In  the  present  analysis 
of  thick  plates  it  was  found  that  there  was  no  need  to  selectively  integrate  the  shear  terms 
over  the  bending  terms.  The  final  analysis  discussed  in  the  results  section  was  done  with 
a  two  point  integration  on  the  bending  terms  and  a  three  point  integration  on  both  the 
mass  and  shear  stiffness  terms. 

The  last  step  in  developing  the  final  form  of  the  global  stiffness  and  mass  matrices  is 
the  elimination  of  the  fixed  degrees  of  freedom  to  implement  the  boundary  conditions.  For 
each  fixed  degree  of  freedom,  the  corresponding  row  and  column  are  eliminated  from  [K] 
and  [M\.  With  this  performed,  eqn  (  2.59)  is  ready  to  be  solved. 

2.4.2  Solution  to  the  eigenvalue  problem 

In  the  last  section  we  discussed  the  integration  of  eqn  (  2.57)  and  eqn  (  2.58)  to  establish 
the  global  matrix  equation,  (eqn  (  2.59)).  The  eigenvalues  and  eigenvectors  of  this  equation 
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provide  the  natural  frequencies  and  mode  shapes  of  the  modeled  plate.  They  are  found 
computationally  using  an  appropriate  linear  algebra  technique.  For  this  task,  the  present 
investigation  utilized  the  IMSL-10  scientific  package  which  is  readily  available  on  the 
CYBER  990  computer.  The  IMSL-10  package  has  a  routine  which  solves  the  general 
eigenvalue  problem  for  the  case  when  [Af]  is  positive  definite,  as  we  have  for  this  problem. 
The  routine  is  efficient  and  returns  till  of  the  desired  eigenvalues  and  their  corresponding 
eigenvectors.  The  eigenvalues  are  ordered  from  lowest  to  highest,  but,  if  under  integration 
is  being  used  for  the  transverse  shear  terms,  care  must  be  taken  in  choosing  the  first 
fundamental  frequency.  In  investigating  the  effects  of  under  integrating  these  terms,  one 
finds  that  it  causes  both  kinematic  and  zero  energy  deformation  modes  to  be  present  in 
the  solution.  Thus,  all  zero  (or  near  zero)  frequencies  need  to  be  overlooked,  as  well  as  any 
whose  deformation  carries  a  zero  energy  deformation  mode  in  any  of  its  elements.  The 
zero  energy  deformation  modes  can  easily  be  spotted  by  looking  at  the  mode  shapes  for 
the  frequency  in  question. 

2.4.3  Updated  natural  frequency  calculation 

As  stated  in  Section  1.6,  the  purpose  of  this  study  is  to  develop  a  technique  to  accurately 
calculate  the  shear  correction  coefficients  so  that  we  c;  perform  an  accurate  vibration 
analysis  of  a  composite  structure.  Having  accurate  values  of  the  shear  correction  coeffi¬ 
cients  allows  us  to  account  for  the  appropriate  transverse  shear  contribution  in  our  strain 
energy  expressions,  and  thus  we  can  receive  accurate  results.  The  first  order  shear  defor¬ 
mation  theory  establishes  a  constant  value  of  transverse  shear  strain  through  the  thickness 
of  the  plate.  In  reality,  it  is  well  known  that  this  is  not  the  case.  In  an  isotropic,  rect¬ 
angular  cross-section  beam  or  plate,  the  exact  elasticity  solution  establishes  a  parabolic 
distribution.  These  two  strain  distributions  result  in  different  transverse  strain  energy 
density  distributions,  which  can  then  be  used  to  calculate  the  shear  correction  factors. 
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This  is  exactly  the  technique  which  will  be  used  in  this  analysis,  not  through  an  exact 
method,  but  by  an  approximate  numerical  method.  We  will  calculate  the  strain  energy  per 
unit  area  of  the  plate  based  on  the  first  order  theory  and  compare  it  to  the  strain  energy 
calculated  using  another  more  accurate  technique.  Both  methods  will  use  the  preliminary 
mode  shape  and  natural  frequency  obtained  using  the  finite  element  analysis.  The  values 
received  for  the  shear  correction  coefficients  should  primarily  be  functions  only  of  the  plate 
geometry  and  lay-up.  It  should  not  be  a  function  of  the  deformation  or  location  in  the 
plate.  The  coefficients  should  be  a  constant  for  a  given  plate,  but  kx  cam  be  different  from 
ky  due  to  inherent  differences  in  the  material  properties  for  the  two  directions. 

The  transverse  strain  energy  per  unit  area  can  be  accurately  calculated  if  one  knows  the 
transverse  stress  distributions  through  the  thickness  of  the  plate,  along  with  the  plate’s 
constitutive  properties.  In  this  analysis  the  transverse  stress  distribution  through  the 
thickness  of  the  plate  will  be  calculated  by  integrating  the  equations  of  equilibrium  with 
respect  to  the  out-of-plane  direction.  Since  the  in-plane  stresses  are  accurately  calculated 
using  the  first  order  analysis,  we  can  expect  to  get  reasonable  results  for  the  transverse 
stresses.  Thus,  with  this  approach  we  can  find  accurate  values  for  the  shear  correction 
coefficients  and  ultimately  the  natural  frequencies  and  mode  shapes. 

2.4.4  The  equations  of  equilibrium 

The  equations  of  equilibrium  for  three  dimensional  elasticity  are  given  by: 


4"  & xy,y  "k  &xz,z  "k  fx 

=  pu 

"k  ayy,y  "k  &yz,z  +  fy 

=  pv 

(2.60) 

+  azy,y  +  C?zz,z  +  fz 

=  pw 

where  /;  are  any  body  forces  present.  The  first  two  terms,  along  with  the  fourth  term  on 
the  left  hand  sides  can  be  brought  over  to  the  right  and  these  forms  then  integrated  with 
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respect  to  z.  The  results  are: 


<Jx  z{z) 


<Tyz(z) 


&zz(z) 


fx)  dz 


fy) dz 

fz)  dz 


(2.61) 


The  first  two  of  these  equations  provide  us  with  a  method  to  calculate  the  transverse 
stresses  axz  and  ayz.  Once  this  is  done,  the  third  provides  us  with  an  expression  for  ozz, 
if  so  desired.  Obviously  to  perform  the  integration  through  the  thickness  of  the  plate  we 
must  first  find  the  in-plane  stresses  and  then  their  derivatives. 

2.4.5  In-plane  stress  calculations 


Once  we  have  found  the  first  eigenvalue  and  eigenvector  of  eqn  (  2.59),  we  know  all  of 
the  nodal  displacements  of  the  plate  model  within  a  constant.  This  constant  will  prove 
to  be  inconsequential  later,  as  we  will  be  interested  only  in  ratios  of  the  strain  energy 
densities.  The  stresses  within  the  plate  can  be  found  from  the  displacements  through  eqns 
(  1.7)-(  1.12)  and  eqn  (  2.4).  Once  these  are  known  at  various  locations,  the  derivatives 
follow. 

From  a  computational  point  of  view,  the  stresses  can  be  calculated  at  a  point  within  an 
element  based  upon  the  element’s  nodal  displacements.  The  usual  procedure  is  to  perform 
the  differentiation  of  the  displacements  by  differentiating  the  shape  function  representation 
of  the  displacements  and  scaling  them  to  the  global  coordinates  through  the  Jacobian. 
The  location  of  the  most  accurate  stress  calculations  can  be  shown  to  be  at  the  Gauss 
integration  points.  Once  the  stresses  are  known  at  the  Gauss  points,  their  derivatives 
can  be  found  in  a  similar  manner.  In  Section  2.3.1  we  briefly  discussed  a  drawback  to 
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the  quadratic  element  chosen  for  this  analysis,  and  it  is  now  clear  that  the  problem  will 
be  in  calculating  accurate  derivatives  of  the  stress  within  an  element.  The  displacements 
within  an  element  will,  at  best,  be  quadratic  and  after  differentiation  will  provide  us  with 
linear  stresses  (in  certain  directions).  The  required  additional  differentiation  leaves  us  with 
constant  stress  derivatives,  and  hence,  constant  strain  energy  per  unit  area  throughout 
the  element  (again  in  certain  directions).  As  a  result,  we  will  have  a  mismatch  when  we 
compare  this  strain  energy  to  the  strain  energy  calculated  using  the  first  order  theory, 
which  results  in  a  higher  order  function.  To  alleviate  this  problem,  a  technique  has  been 
developed  to  globally  smooth  the  nodal  displacements  to  a  higher  order  polynomial.  When 
the  new  displacements  are  used  in  conjunction  with  higher  order  shape  functions,  the  end 
result  is  a  higher  order  variation  of  transverse  strain  energy  within  the  element.  We  will 
find  that  this  linear  variation  will  be  sufficient  in  calculating  the  new  shear  correction 
coefficients. 

The  method  is  begun  by  considering  each  element  individually.  Using  a  cubic  least 
squares  curve  fit  routine,  we  calculate  the  four  polynomial  coefficients  for  w,  <f>x  and  <f>y  on 
each  side  of  the  element.  In  other  words,  we  end  up  with  the  coefficient  data  for  twelve 
curves  for  each  element.  The  curve  fitting  is  done  utilizing  five  points  per  side.  We  use  the 
original  three  data  points  plus  two  additional  ones.  If  the  element  is  an  interior  element, 
a  data  point  from  either  side  is  taken.  If  it  is  a  boundary  element,  two  points  from  the 
adjacent  interior  element  jure  used.  After  this  is  completed,  the  next  step  is  to  temporarily 
transform  our  eight  noded  element  into  a  twelve  noded  one,  as  depicted  in  Figure  2.3.  The 
shape  functions  associated  with  this  element  are  established  as: 
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Figure  2.3:  Twelve  Noded  Isoparametric  Element 


Ni  =  &  (1  +  «<)  (1  +  W.)  (9^2  +  V  -  10)  *  =  1, 2, 3, 4 

Ni  =  £  (1  +  ££<)  (1  +  9777/,)  (1  -  t?2)  t  =  7,8, 11, 12  (2.62) 

iV,  =  ^  (1  +  9ft.)  (1  +  Wl)  (1  -  £2)  t  =  5, 6, 9, 10 

Next,  we  proceed  to  build  an  element  displacement  vector,  {A}c,  consisting  of  sixty  terms 
and  defined  as 


{A}c  = 


(«o)l 

Kh 

K)i 

(<Pz)l 

(4>y)l 


[  (4>y)  12  J  c 


(2.63) 


Note  that  the  c  subscript  distinguishes  this  displacement  vector  from  the  original  one 
defined  in  eqn  (  2.32).  In  building  {A}c,  the  original  displacements  axe  retained  for  the 
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comer  nodes,  and  the  values  for  the  side  nodes  are  calculated  using  the  polynomial  curve 
fits.  In  this  manner,  the  displacements  which  originally  varied  quadratically,  now  vary 
cubically. 

Next,  the  stresses  are  calculated  at  the  twelve  nodes.  We  start  by  finding  the  strains 
from  the  derivatives  of  displacements  at  each  node  by  establishing 

{ £; }  «  m  {a}.  (2.64) 


where  [/?*]  is  defined  exactly  as  [/3]  was  defined  in  eqns  (  2.28)-(  2.35)  with  the  exception 
that  [III]  is  expanded  out  to  include  the  extra  four  shape  functions.  It  is  important  to 
remember  that  eqn  (  2.64)  provides  the  midplane  strains  and  curvatures  at  a  specific  point 
in  the  element.  The  three  in-plane  strains  at  this  point  can  be  represented  as  a  function 
of  z  by 
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The  stresses  at  this  specific  point  in  the  plate’s  thickness  become 


(7  xz  | 

’  °yy  f 
,  °*v  ) 


[Q]k  [Z]  r ]  {A}c 


(2.65) 


(2.66) 


where  [Q\k  represents  the  [Q],  as  defined  in  eqn  (  2.9),  for  the  fcth  layer  of  the  laminate, 
and  z  falls  in  this  layer.  Thus,  eqn  (  2.66)  allows  us  to  calculate  the  stress  at  a  specific 
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r],  £,  z  location  in  the  element.  This  is  the  first  necessary  step  in  calculating  the  derivatives 
of  the  stresses  in  the  element. 

2.4.6  Derivatives  of  stresses 


In  Section  2.4.4  we  established  that  in  order  to  calculate  the  transverse  stresses,  and 
ultimately  the  transverse  strain  energy  in  the  laminated  plate,  we  needed  to  calculate 
the  derivatives  of  the  in-plane  stresses  with  respect  to  the  in-plane  directions.  Above,  in 
Section  2.4.5,  we  established  a  method  to  calculate  the  in-plane  stresses  at  any  location 
within  the  element  using  eqn  (  2.66).  This  equation  will  become  the  basis  for  the  next 
calculations. 

We  start  by  choosing  a  convenient  set  of  points  at  which  to  determine  our  desired 
parameters.  Through  experimentation  and  ease  of  computational  implementation,  it  was 
found  best  to  use  eight  interior  points  coinciding  to  the  Gauss  integration  points  for  a 
three  by  three  integration,  (the  center  point  is  discarded).  The  stresses,  as  given  by  eqn 
(  2.66),  are  calculated  as  a  function  of  z  at  each  of  these  eight  interior  points.  The  result 
is  established  in  matrix  notation  as: 


(  (")i  ’ 
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(  (*)8  ) 
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(2.67) 


Turning  our  attention  now  to  the  derivatives  of  the  stresses,  we  establish  only  the  terms 


needed  and  write  them  as 
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(2.68) 


where 


ru  r  12 
T21  T22 


(2.69) 


The  column  matrix  on  the  right  hand  side  of  eqn  (  2.68)  can 


now  be  established  through 


eqn  (  2.67)  and  the  derivatives  of  eqn  (  2.62).  The  result  is 
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(2.70) 


We  can  combine  eqn  (  2.68)  and  eqn  (  2.70)  to  get  an  expression  for  the  derivatives  of 
stresses  at  any  rj,£  as  a  function  of  z.  The  result  is: 
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(2.71) 


where  the  rectangular  matrices  in  eqn  (  2.68)  and  eqn  (  2.70)  have  been  represented  by 
[G]  and  [AT]  respectively. 
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Body  force  or  dynamic  terms 

The  representation  of  the  body  forces,  or  the  dynamic  terms,  is  the  final  item  required 
before  the  integration  of  eqn  (  2.61)  can  be  performed.  For  the  case  of  free  vibration, 
the  inertia  terms  of  the  plate  can  be  treated  either  way  with  no  difference,  as  the  sign 
difference  drops  out.  Recalling  the  discussion  in  Section  2.3.3  we  can  express  the  required 
acceleration  terms  as: 

u  —  —u2  ( uQ  +  z<f>x) 

v  =  — cj2  ( va  +  z<t>y)  (2.72) 

2 

w  =  —uw0 

These  equations  are  the  final  expressions  we  need  to  calculate  the  transverse  stress  distri¬ 
butions. 

Putting  it  all  together 

Equation  2.71  and  eqn  (  2.72)  provide  the  appropriate  expressions  to  place  into  eqn  (  2.61) 
so  that  the  integration  can  be  accomplished.  The  integration  can  be  performed  numerically 
using  a  trapezoidal  rule,  or  other  simple  integration  method,  starting  at  z  =  —  h/2  with 
azz  and  ayz  equal  to  zero,  or  to  whatever  surface  tractions  exist.  The  final  values  of  axz 
and  (Tyz  at  z  =  h/2  should  correspond  to  the  tractions  on  the  top  surface  of  the  plate, 
which  are  zero  for  the  free  vibration  case.  The  end  results  of  the  calculations  are  the 
determinations  of  the  transverse  bhear  distributions  through  the  thickness  of  the  plate  at 
some  specific  location. 


48 


2.4.7  Shear  correction  coefficient  calculations 
First  order  strain  energy  calculations 

The  transverse  strain  energy  densities  calculated  using  the  first  order  theory  are  Uxz  and 
Uyz.  They  are  defined  by 


Uxz  =  ~Cxzyl 
Uyz  ~  2  U'yz'^yz 

where 


(2.73) 


(2.74) 


and  yxz  and  yyz  are  given  in  eqns  (  1.10)-(  1.11). 


Improved  strain  energy  calculations 

As  discussed  earlier,  through  the  use  of  the  accurate  in-plane  stresses  and  the  equilibrium 
equations,  we  can  calculate  a  more  accurate  representation  of  the  transverse  strain  ener¬ 
gies.  The  transverse  stress  distributions,  as  functions  of  z,  are  found  by  substituting  eqn 
(  2.71)  and  eqn  (  2.72)  into  eqn  (  2.61).  The  result  of  this  provides  us  with  through-the- 
thickness  distributions  of  both  axz  and  oyz.  The  transverse  strain  energies  due  to  these 


stresses  are  defined  as: 
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1  rh/ 2  _  /  \ 

U„  =  -  /  C44  hi)  dz  (2.75) 

2  J-/1/2  v  ' 

1  rh /2  _  /  „  \ 

^  =  2].knC^-)dz  (2'76) 

where 


Tarr  —  ^44^xz  "h  S^ffyz 

T y2  =  ^54^zz  "h  ^55^yz 


Here 


S44  ^45 

Q44  Q 45 

S54  S55 

Q54  Q55 

(2.77) 


Transverse  strain  energy  comparisons 


In  the  previous  two  sections  we  have  developed  expressions  for  the  transverse  strain  energy 
using  the  results  from  a  first  order  theory  (eqns  (  2.73)-(  2.74))  and  by  integrating  the 
equations  of  equilibrium  through  the  thickness  of  the  plate  (eqns  (  2.75)-(  2.76)).  By 
comparing  the  results  of  the  two  different  methods  we  can  obtain  a  new  shear  correction 
factor  with  which  an  updated  natural  frequency  can  be  calculated.  If  we  let  k°  and  k°  be 
the  original  shear  correction  coefficients,  then 


(k,)2Uiz  =  ( k°)2Ulz 


(2.78) 


should  hold  true,  and  hence  new  shear  correction  coefficients  can  be  found. 


CHAPTER  III 


LEAST  SQUARES  ELEMENT 
DEVELOPMENT:  Cl  CONTINUITY 
APPROXIMATION  THROUGH  A  LEAST 
SQUARES  METHOD 

3.1  Background  Information 

As  discussed  in  Section  1.3,  one  approach  to  improving  laminated  plate  theory  is  to  assume 
a  displacement  field  with  higher  order  terms  in  the  z  coordinate.  For  example,  a  typical 
form  often  seen  (see  eqns  (  1.15))  assumes  «  and  v  cubic  in  z  and  w  quadratic.  Rewriting 
the  displacement  field  here  for  convenience,  we  have: 

u  =  ua  +  z<j>x  +  z2ipx  +  z3(x 

v  =  v0  +  z<f)y  +  z2il>y  +  z3Cy  (3.1) 

w  =  w0  +  Z(f>z  +  z2ipz 

The  motivation  behind  cubic  functions  for  u  and  v  stems  from  the  fact  that  transverse 
shear  should  at  least  be  quadratic  in  z.  This  allows  the  top  and  bottom  surfaces  to  be 
free  from  tractions,  if  conditions  warrant  this.  The  expression  for  w  is  generally  simplified 
by  dropping  the  terms  involving  z,  thus  making  w  a  constant.  This  displacement  field  has 
been  shown  to  be  an  improvement  over  the  first  order  theory,  mainly  for  the  reason  of 
producing  acceptable  results  without  the  use  of  shear  correction  factors.  The  disadvantage 
is  that  the  displacement  field  is  now  an  eleven  parameter  field  as  shown  above,  or  assuming 
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w  as  a  constant,  a  nine  parameter  field.  When  compared  to  the  five  parameter  first 
order  theory,  this  difference  represents  a  significant  increase  in  the  number  of  unknown 
functions  in  the  problem.  The  nine  parameter  field  can  easily  be  simplified,  as  has  been 
done  beginning  in  1984  by  Reddy  [113],  by  setting  the  transverse  stress  to  zero  on  the  top 
and  bottom  surfaces  of  the  plate.  These  four  conditions  are  used  to  eliminate  the  variables 
V’z>Vv>Cr  and  Cy  from  the  u  and  v  displacements.  The  end  result  is: 


u  =  u„  +  z 


4z2 

4>x  ~  x  d-  wo,x) 


v  =  vQ  +  z 


4  z2 

4>y  ~  (fiy  d"  wo,y) 


(3.2) 


w  =  w0 

This  field  is  readily  seen  to  now  have  five  unknowns,  which  is  the  same  number  as  the 
first  order  theory.  This  displacement  field  provides  very  good  results  when  applied  to 
laminated  plates  and  is  considered  an  improvement  over  the  first  order  theory.  It  performs 
well  without  shear  correction  coefficients  and  provides  improved  in-plane  responses,  and 
hence,  can  be  used  to  provide  better  out-of-plane  results  following  the  method  outlined  in 
Section  2.4.7.  This  improvement  does  not  come  without  cost,  however.  Even  though  the 
number  of  unknown  functions  is  the  same  as  in  the  first  order  theory,  the  finite  element 
formulation  of  this  field  is  hampered  by  the  form  of  eqns  (  3.2),  which  include  derivatives 
of  w  in  the  u  and  v  displacements.  In  other  words,  the  terms  w,z  and  w.y,  found  in 
the  expressions  for  u  and  v,  make  the  problem  not  one  of  C°  continuity,  as  it  was  for 
the  first  order  and  higher  order  theories,  but  make  it  now  one  requiring  C1  continuity 
for  w.  This  requirement  is  the  reason  that  the  displacement  field  in  eqns  (  3.2)  is  not 
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widely  used.  The  C1  continuity  problem  has  long  plagued  finite  element  developers  and 
is  usually  avoided  if  at  all  possible.  In  fact,  one  of  the  biggest  advantages  of  the  five 
parameter  Mindlin  style  displacement  field  was  that  it  changed  the  plate  bending  problem 
from  one  of  C1  continuity  found  in  the  Kirchhoff  theory,  to  one  of  C°  continuity.  The 
standard  isoparametric  shape  functions,  given  in  eqn  (  2.21),  insure  that  the  degrees  of 
freedom  are  continuous  from  element  to  element  but  do  nothing  to  maintain  continuity 
in  the  derivatives  of  the  functions.  In  most  analysis  this  fact  only  becomes  significant 
when  differentials  of  the  degrees  of  freedom  are  used,  for  instance,  as  required  in  stress 
calculations.  The  stresses  are  discontinuous  between  elements  because  of  the  discontinuity 
in  the  shape  function  differentials  from  element  to  element.  In  the  case  of  eqns  (  3.8), 
the  discontinuity  of  w,x  and  w,y  immediately  cause  discontinuities  in  u  and  v.  Thus,  the 
global  displacement  field  contains  gaps,  violating  one  of  the  requirements  for  convergence1. 
The  result  is  a  formulation  which  is  generally  too  soft  and  has  poor  convergence  qualities. 
In  fact,  a  displacement  field  of  the  form  given  in  eqn  (  3.2)  was  implemented  into  a  C° 
element  formulation  during  this  course  of  study  with  very  poor  results. 

Successful  use  of  eqns  (  3.2)  in  a  finite  element  formulation  has  been  demonstrated  by 
Reddy  and  Phan  (1985)  [105],  but  formulations  of  this  type  are  by  no  means  common. 
The  problem  was  solved  by  utilizing  Hermite  polynomials  for  the  element  shape  functions 
which  satisfy  the  C1  continuity  requirement.  However,  the  practicability  of  Hermitian 
polynomials  is  questionable  based  upon  the  lack  of  published  work  utilizing  them.  In 
fact,  most  finite  element  text  books  have  little  or  no  reference  to  them.  See  Zienkiewicz 
[146],  Bathe  [9]  and  Cook  [23].  In  the  book  by  Zienkiewicz,  he  states  that  elements  using 
Hermitian  interpolation  functions  have  little  engineering  applicability. 

In  the  past,  researchers  have  proposed  several  other  methods  to  solve  the  C1  continuity 
1  It  can  be  shown  that  convergence  can  still  be  achieved  if  gaps  go  to  zero  in  the  limit 
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problem.  One  common  method  discussed  by  Zienkiewicz  has  been  the  use  of  non-conformal 
elements,  which  allow  gaps,  but  satisfy  convergence.  This  method  has  been  used  with  some 
success  in  other  problems,  but  it  immediately  restricts  the  elements  to  being  straight  sided 
rectangles.  Another  possible  solution  includes  the  use  of  conforming  triangular  elements  as 
developed  by  Anderheggen  (1970)  [4]  and  Irons  (1969)  [37].  These  have  been  demonstrated 
for  thin  plates. 

All  of  the  above  mentioned  techniques  to  solve  the  C 1  continuity  problem  generally  are 
not  without  their  drawbacks.  It  is  the  intent  of  this  research  to  present  a  simple,  unique 
and  computationally  effective  method  to  approximate  a  C1  continuous  element  which 
produces  accurate  results  for  laminated  composite  plates.  In  the  next  several  sections,  the 
new  Least  Squares  element  will  be  developed  in  detail  and  will  approximate  C1  continuity 
with  very  little  increase  in  complexity  over  traditional  finite  element  techniques.  We  will 
show  that  fields  of  the  form  of  eqns  (3.2)  can  be  used  easily  with  very  good  results. 

3.2  Displacement  Field  Development 

3.2.1  Displacement  Field  Basis 

For  this  work  we  will  develop  a  displacement  field  similar  to  that  given  in  eqns  (  3.2), 
but  which  will  be  better  suited  for  our  needs  later  in  developing  a  piecewise  continuous 
displacement  field.  We  begin  by  assuming  a  symmetric,  parabolic,  transverse  strain  field 
of  the  form: 


laz  =  <Pa(ao  +  a\z 


2), 


(3.3) 


where  is  a  transverse  shear  strain,  a  =  x.y  and  ao  and  ai  are  constants.  If  we  force 
the  transverse  stress  (strain)  to  go  to  zero  at  z  =  ±^,  we  can  eliminate  the  two  constants. 
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The  result  can  be  written  as: 


Next,  substituting  eqn  (  3.4)  into  the  strain  displacement  relations  for  transverse  strain, 
(see  eqns  (  1.10)-(  1.11)),  and  solving  for  uQ)Z,  we  can  write: 


Integrating  this  with  respect  to  z  yields 


«Q  =  «QO  +  VaU-^2 


-  ZW,Q 


The  end  result  is  a  displacement  field  of  the  form: 


U  =  Uo  +  VxyZ-y-p]-  zu>,x 


(  4*3 ) 

v  =  vo  +  <Py  U-  3^2  )  ~ZW'V 


w  =  w(x,y) 


A  displacement  field  of  this  form  was  published  by  Bhimaraddi  and  Stevens  [16]  in  1984, 
but  it  was  presented  with  no  justification  or  explanation  of  how  or  why  it  was  chosen. 
The  displacement  field  is  somewhat  similar  to  that  given  in  eqns  (  3.2),  but  instead  of  a 
rotation  angle,  4>a,  we  have  a  measure  of  transverse  shear  strain, 
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The  field  in  eqns  (  3.7)  has  no  immediate  advantage  over  eqns  (  3.2),  as  they  both 
require  C1  continuity.  In  fact,  eqns  (  3.7)  are  at  somewhat  of  a  disadvantage  over  eqns 
(  3.2)  if  one  intends  to  compare  any  results  to  those  of  a  first  order  shear  deformation 
analysis  because  of  the  differences  between  <j>a  and  <pQ.  The  advantage  to  eqns  (  3.7)  comes 
into  play  when  it  is  used  as  a  basis  for  a  simplified  piecewise  continuous  displacement 
model.  The  algebra  involved  in  satisfying  continuity  becomes  much  simpler.  Therefore, 
throughout  this  section,  eqns  (  3.7)  will  be  utilized  as  the  displacement  field  of  choice. 

3.2.2  The  Piecewise  Continuous  Displacement  Field 

As  discussed  earlier  in  Section  1.4.2,  a  continuous,  smooth  function  cannot  accurately 
represent  the  displacements  through  the  thickness  of  a  laminated  composite  plate.  To 
improve  the  accuracy  of  the  analysis,  we  must  turn  to  a  piecewise  continuous  function 
through  the  thickness  where  each  layer  is  modeled  with  its  own  smooth  function.  This 
allows  for  slope  discontinuities  at  the  layer  interfaces.  Towards  this  end,  we  assume  a 
function  of  the  form  given  in  eqns  (  3.7),  but  to  begin  with,  we  allow  each  layer  to  have 
its  own  values  for  u0,v0,<px  and  <py.  This  cam  be  written  as: 


wn  =  w(x,y) 

where  the  superscript  n  refers  to  the  nth  layer  of  the  laminate,  and  2  is  a  global  coordinate 
running  from  —  |  to  4-|. 

The  piecewise  displacement  field  presented  in  eqns  (  3.8)  must  be  modified  to  meet 
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some  of  the  requirements  listed  in  Section  1.4.2,  namely  transverse  traction  continuity  and 
displacement  continuity  at  the  layer  interfaces.  Using  eqns  (  3.8),  we  can  calculate  the 
transverse  stresses  on  either  side  of  the  layer  interfaces  and  set  them  equal  to  one  another. 
In  doing  so  we  find  out  that  the  values  for  can  all  be  related  to  the  ipa  of  one  reference 
layer.  If  we  choose  the  bottom  layer  as  the  reference  layer,  we  can  wnte: 

*>;  =  el  =  ',,eel  (3.8) 

*>;  =  (3.io) 

where  Cl \  and  C55  are  from  eqn  (2.5)  for  the  nth  layer.  Similarly,  enforcing  u  and  v  to 
be  continuous  functions  across  the  layer  interfaces  results  in  the  expressions 

K  =  «o  +  ^2  (Zk  ~ 

k= 1  \ 

<  =  »!  +  EU- 

k=l  \ 

where  the  bottom  layer  is  once  again  the  reference  layer.  With  these  expressions,  eqns 
(  3.8)  can  be  written  as: 


|§V-. -&)*>; 


(3.12) 


4 zjt  \  x 

^ 2  )  -aC<Px 


(3.11) 


u 


n 


K  + 


<p\  +  an<plx  ^ 


4  *3\ 

Mj)-ZW 


v 


n 


+ 


4  z*\ 

3V)-ZW 


(3.13) 
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wn  =  w(x,y) 

v.  here 

Pj  =  (zj  -  4z*/2h2}  {a j-i  -  otj) 

Pj  =  (zj-Az?/3h2)(0j-l-pj) 

and  here  Zj  is  the  distance  to  the  bottom  of  jth  layer.  This  piecewise  continuous,  simplified 
higher  order  displacement  field  will  be  the  one  used  to  demonstrate  the  effectiveness  of 
the  Least  Squares  element  in  the  sections  to  follow.  A  field  of  this  form  was  used  by  Lee 
et  al  (1990)  [59], 

At  this  time  an  important  distinction  is  made  in  terminology.  The  three  equations 
above  will  be  referred  to  as  the  displacement  fields,  while  the  five  terms  on  the  right  hand 
side,  ua,  v0,  ipx,  <py  and  w,  will  be  termed  the  displacement  functions.  This  terminology 
will  help  make  the  following  sections  less  confusing.  Note  also  that  the  superscript  1  has 
been  dropped.  It  will  be  implied  from  this  point  on. 

3.3  Theory  Development:  Method  I 

3.3.1  The  Domain  Displacement  Functions 

We  begin  by  assuming  the  element  domain  shown  in  Figure  3.1  with  a  local  r?-£  coordinate 
system  as  shown.  Within  this  domain  we  have  displacements  represented  by  eqns  (  3.13)  in 
terms  of  five  domain  displacement  functions  u0,  va,  ui,  <pz  and  (py.  Each  of  these  functions 
is  represented  by  an  n  term  polynomial  expansions  in  rj  and  £  with  n  unknown  constants 
or,.  The  expressions  can  be  written  as: 

u0  =  £*[  +  0127]  +  a3£  +  a^T]  +  +  a ^  - 
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Figure  3.1:  Least  Squares  Element  Domain 

v0  =  an+\  +  a„+2 7?  +  an+3^  +  «n+4^?2  +  <*n+5s2  +  an+6^  H - 

V^X  =  «2n+l  +  OC2n+2rJ  +  «2n+3^  +  «2r»+47?2  +  <*2n+5£2  +  a2n+6Vt  H (3-14) 

V’tf  =  «3n+l  +  «3n+2^?  +  <*3r»+3£  +  C*3n+4f?2  +  a3n+5^2  +  UZn+eVt  4 - 

w  =  ®4n  +  l  +  tt4n+27?  +  a4n+3£  +  a4r,+47/2  +  a4n+5£2  +  a4n+6T/£  H - 


We  can  write  these  expressions  in  matrix  notation,  (using  n  =  10  for  example),  as: 


«o 

.  «10 

Vo 

:  : 

<Px 

*  = 

<Py 

*  *  *  * 

.  W  . 

«41  .  a, 50 

{A} 


(3.15) 


where 
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{A}t  =  [AJ  =  LI  T]  Z  T]Z  r?  Z2  V2t.  vt2  v3  f3J 


Or  we  can  write  the  functions  in  terms  of  a’s  as: 


where 


Vo  | 

-  <Px  ► 
<fiy 


W  M 


Ml 


'AO  0  0  O' 

0  A  0  0  0 

0  0  A  0  0 

0  0  0  A  0 

.  0  0  0  0  A 


and  {a}  is  a  column  vector  of  all  the  a’s  (50  for  this  case).  Next,  we  add 
derivatives  of  w  to  these  expressions  to  get: 


where 


u0 

v0 

<Pz 

'  Vy  • 
w 
IV,  T, 

.  W<  . 


[A]  {«} 


Ml 


A  0  0  0  0 

0  A  0  0  0 

0  0  A  0  0 

0  0  0  A  0 

0  0  0  0  A 

0  0  0  0  A.v 

0  0  0  0  A,€ 


(3.16) 


(3.17) 


(3.18) 


the  local 


(3.19) 


(3.20) 
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We  now  have  everything  we  need  to  write  the  domain  displacement  field  given  in  eqns 
(  3.13)  in  terms  of  the  domain  functions  given  in  eqns  (  3.14).  The  result  is 


Here 


=  PI  [A]  {a} 


where 


1  0  ci  0  0  C2Tii  c2ri2 

o  i  o  ci  o  c2r2i  c2r22 

0  0  0  0  1  0  0 


(3.21) 


(3.22) 


ci 


n 


Epi 


+  zatk  - 


4z3 

3&ak 


(3.23) 


ci 


(3.24) 


c2 


4  z3 
3  h2 


(3.25) 


This  now  gives  an  expression  for  the  domain  displacement  field  in  terms  of  the  a’s  .  The 
next  step  is  to  develop  a  boundary  displacement  field  written  in  terms  of  nodal  degrees  of 
freedom. 


3.3.2  Boundary  Displacement  Functions  from  Nodal  Degrees  of  Free¬ 
dom 


The  element  domain  shown  in  Figure  3.1  is  now  modified  to  include  eight  boundary  nodes. 
The  element  is  shown  in  Figure  3.2  and  now  resembles  a  standard  quadratic  isoparametric 
element.  However,  we  define  the  nodal  degrees  of  freedom  as 
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{6,}T  = 

[  Hq 

V0  Vx  <Py  w 

w,x  w,y  J  i  =  1,2, 3,4 

(3.26) 

{s,}T  = 

[  u0 

Vo  Vx  V>y  J 

i  =  5,6,  7, 8 

(3.27) 

giving  the  side  nodes  and  the  corner  nodes  different  numbers  of  degrees  of  freedom  for  a 
total  of  44  degrees  of  freedom  per  element.  For  future  use,  we  define  a  vector  containing 
all  44  degrees  of  freedom  as: 


(3.28) 


Along  each  edge  of  the  element  we  have  two  comer  nodes  with  7  degrees  of  freedom 
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n 


s  =  -1.0  s  =  0.0  s  =  +1.0 


Figure  3.3:  Variation  of  Nodal  variable  q  along  element  side. 

and  one  side  node  with  only  four  degrees  of  freedom,  as  described  above.  With  this,  uoy 
vQ,  <pz  and  i py  are  defined  by  three  quantities  along  the  edge,  one  at  each  node.  Any  one 
of  these  four  quantities,  call  it  q,  defines  a  quadratic  curve  in  the  boundary  variable  s, 
which  varies  from  —1.0  to  +1.0  along  each  of  the  four  sides.  See  Figure  3.3.  We  can 
easily  establish  the  second  order  equation  for  any  of  the  four  degrees  of  freedom  from  the 
following  equation: 

?(»)  =  £  (~s  +  s2)  ?(_)  +  (l  -  s2)  q(0)  +  ^  (s  +  s2)  q(+)  (3.29) 

The  function  w ,  on  the  other  hand,  is  defined  by  six  nodal  variables.  We  define  first  w,s 
as  the  derivative  of  w  tangent  to  the  side  in  the  direction  of  s,  and  w,n  as  the  derivative 
of  w  normal  to  the  direction  of  s.  See  Figure  3.4.  For  a  rectangular  element  aligned  with 
the  global  x-y  coordinate  system  these  two  derivatives  correspond  identically  to  either  w,z 
or  w,y  depending  upon  which  side  one  looks  at.  If  the  element  is  not  straight  sided  or 
not  aligned  with  the  global  coordinate  system,  the  two  derivatives  will  each  be  a  linear 
combination  of  both  w,z  and  w.y  related  through  the  standard  Jacobian  matrix.  With 
this  in  mind,  we  see  that  w  is  defined  by  four  variables,  the  two  values  of  w  and  the  two 
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Figure  3.4:  Variation  of  Nodal  variable  w  along  element  side. 

values  w,a.  Thus,  w  can  be  fit  by  a  cubic  equation,  making  w,a  quadratic.  In  addition, 
w,n  is  defined  only  by  two  points,  one  io,„  at  each  end,  thus  making  it  vary  only  linearly. 
The  required  equations  for  w  and  its  derivatives  become: 


(s)  =  Q  ~  +  ^s3)  wi  }  +  \  (x  ~  s  ~  s2  +  s3)  ] 

+  Q  +  IS  ~  ^+)  +  \  (-1  -  s  +  s2  +  s3)  «,<+>  (3.30) 


2 

W,s{s)  =  ~  (-1  +  S2) 

|  wf'  )  +  \  1 

4 

(-1  -2s  +  3s2) 

i  ™(+)  +  j  i 

4 

(-l  +  2a  +  3s2) 

«V.M  =  Q  -  |)  +  Q  +  |)  <3-32> 


The  significant  point  here  is  that  all  functions,  including  w„  and  to,„,  can  be  expressed 
along  an  edge  of  an  element  in  terms  of  only  the  nodal  values  along  that  particular  edge. 
This  will  become  an  important  fact  to  help  insure  compatibility  from  element  to  element. 
The  final  step  is  to  compile  cqns  (  3.29)-(  3.32)  into  a  matrix  format.  It  is  convenient 
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to  write: 


u0 

Vo 

<Pz 
•  <Py 
w 
VJ,x 
,  W'V 


=  m,{A} 


(3.33) 


where  the  subscript  i  refers  to  the  ith  side  of  the  element,  and  the  matrix  [!F\  contains 
the  necessary  terms  from  eqns  (  3.29)-(  3.32)  as  well  as  Jacobian  terms  relating  the  local 
differentials  to  the  global  ones. 

3.3.3  The  Least  Squares  Implementation:  Method  I 

In  the  past  two  sections  we  have  developed  two  different  methods  for  finding  expressions 
for  the  five  displacement  functions ,  u0 ,  v0,  <px,  <py  and  w,  as  well  as  for  the  differentials  of 
w.  In  Section  3.3.1,  eqn  (  3.19)  defined  the  functions  anywhere  in  the  element’s  domain 
in  terms  of  unknown  a’s.  In  Section  3.3.2,  eqn  (  3.33)  defined  the  seven  functions  only 
on  the  element’s  boundary  and  in  terms  of  the  nodal  degrees  of  freedom.  Let  us  refer  to 
column  vectors  of  these  seven  functions  from  each  of  these  two  equations  as  {<!>n}  and  {6p} 
respectively.  Here  the  symbols  Cl  and  T  correspond,  of  course,  to  domain  and  boundary. 
The  vector  {6q}  gives  the  seven  functions  in  terms  of  the  local  coordinates,  rj  and  £,  while 
{<5r }  gives  them  as  functions  of  a  boundary  variable,  s. 

The  basis  for  the  formulation  of  the  Least  Squares  element  begins  with  these  two 
vectors.  The  vector  {£p}  is  evaluated  on  the  boundary  by  setting  the  appropriate  variable, 
either  rj  or  £,  to  ±1.0  while  the  other  varies  as  s.  This  new  vector  will  be  referred  to  as 
{5n}r,  in  other  words,  the  domain  displacement  functions  evaluated  on  the  boundary. 
Next,  we  define  a  functional,  /,  as  the  integral  of  the  square  of  the  difference  between  the 
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domain  displacement  functions  evaluated  on  the  boundary  and  the  boundary  displacement 
functions  established  from  the  nodal  degrees  of  freedom.  This  is  written  as: 


/  =  Jr({6n}r-{6r})2dr 


(3.34) 


Substituting  in  the  expressions  from  eqn  (  3.19)  and  eqn  (  3.33)  the  above  equation 
becomes: 


/  =  J^[A]  {<>}-[?],  {A})2  dT 


(3.35) 


The  functional  I  can  now  be  minimized  with  respect  to  the  a’s  to  obtain  the  expression 
1 1-=  f  (Hf  M!  (a)  -  n  (A))  dr  =  0  (3.36) 

(The  constant  2  has  been  divided  out.)  The  next  step  is  to  perform  the  integration  and 
solve  for  the  a’s  in  terms  of  the  nodal  degrees  of  freedom.  As  indicated  in  eqn  (  3.36), 
the  integration  is  performed  over  the  element’s  boundary.  The  integrand  is  defined  as 
a  function  of  the  boundary  variable  s,  but  care  must  be  taken  in  how  the  integration  is 
performed  to  insure  the  correct  sign  on  the  integral  for  each  side.  Upon  performing  the 
integration,  the  expression  can  be  solved  for  the  a’s  2.  The  result  is 

<«>  =  (MTM)~IMTm<A}  (3.37) 

or 

2  Note:  To  keep  the  number  of  symbols  used  to  a  minimum,  the  same  symbols  are  used  after  the 
integration  to  represent  the  variables.  The  difference  is  understood. 
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{«}  =  [HI]{A} 


(3.38) 


where 

{H-\  =  {MT[A\)~'[ATm  (339) 


This  result  is  the  basis  for  the  Least  Squares  Element.  The  unknowns,  (the  ar’s),  of  the 
domain  displacement  functions  which  make  up  the  displacement  fields,  are  chosen  in  such 
a  manner  as  to  force  these  functions  to  match  those  on  the  boundary  which  are  determined 
from  the  nodal  degrees  of  freedom.  Since  the  functions  on  each  side  of  the  element  are 
determined  only  from  the  nodal  degrees  of  freedom  on  that  side,  interelement  compatibility 
is  met.  Note,  in  the  last  two  sentences  the  words  ‘force’  and  ‘met’  were  emphasized. 
This  is  because  the  desired  action  is  only  accomplished  in  an  approximate  manner.  The 
values  of  the  functions  may  not  match  exactly,  but  their  differences  are  minimized.  Thus, 
compatibility  is  not  met  unconditionally,  but  it  is  met  in  a  least  squared  sense. 

With  eqn  (  3.38)  we  now  have  the  ability  to  update  eqn  (  3.19).  Substituting  in  the 
expression  for  {a},  the  domain  displacement  functions  become: 


Va 


V* 

>  <fy 

w 


w"i 

l  w,(  j 


(3.40) 


Remembering  that  these  domain  displacement  functions  form  the  basis  for  the  element 
displacement  field  throughout  the  domain,  we  can  also  update  eqn  (  3.21)  to  now  be: 
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I  >  l  =  [I]M][ff,|{A}  (3.41) 


With  eqn  (  3.41)  we  now  have  an  expression  for  the  element  displacement  field.  This 
displacement  field  formulation  is  quite  different  from  that  of  a  standard  isoparametric 
formulation  and  has  some  points  worthy  of  discussion.  First  of  all,  the  matrix  [*4]  is  a 
function  of  the  elements  local  coordinate  system  variables  r?  and  f .  The  order  of  the  terms 
is  dependant  upon  the  value  of  n  chosen  in  eqns  (  3.14).  Thus,  eqn  (  3.41)  represents 
the  displacement  field  already  in  terms  of  the  element  local  coordinate  system.  The  role 
usually  performed  by  shape  functions  has  already  been  filled  directly  for  the  Least  Squares 
element.  First,  the  strain  field,  or  any  desired  differential  field  for  that  matter,  is  found  by 
directly  differentiating  [A]  in  eqn  (  3.41).  Next,  the  development  of  [Hi]  in  eqns  (  3.35)- 
(  3.38)  has  insured  that  the  domain  displacement  field  has  been  determined  such  that  it 
matches  a  specific  function  determined  from  nodal  degrees  of  freedom  on  each  specific 
edge  of  the  element.  This  specific  function  is  quadratic  for  uQ,  va,  <px  and  <py.  This  is 
essentially  the  same  as  what  is  provided  by  the  shape  functions  in  eqn  (  2.21).  In  fact,  if 
n  in  eqn  (  3.14)  is  chosen  to  provide  a  biquadratic  function  (with  the  t?2£2  term  removed) 
the  formulation  for  these  four  domain  functions  should  be  identical  with  that  found  using 
the  standard  shape  function  interpolations.  However,  unlike  eqn  (  2.21)  ,  the  formulation 
in  eqn  (  3.41)  allows  for  a  cubic  variation  of  w  along  its  edges,  and  hence  a  quadratic 
tangential  derivative,  while  establishing  a  linear  normal  derivative  of  w  along  each  edge. 
This  is  where  the  power  of  the  Least  Squares  element  comes  into  play.  All  three  variables, 
w,  w,x  and  w,y,  will  remain  compatible  along  common  element  boundaries.  As  a  result, 
so  will  u  and  v. 


3.4  Theory  Development:  Method  II 

3.4.1  The  Domain  Displacement  Fields 
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In  Section  3.3.1  we  defined  the  domain  displacement  functions  for  the  element  domain. 
The  unknowns  in  these  expressions  were  then  chosen  to  minimize  the  difference  between 
them  and  the  boundary  displacement  functions  in  terms  of  the  nodal  degrees  of  freedom. 
Method  II  will  utilize  the  displacement  fields  themselves  in  the  Least  Squares  method 
rather  chan  the  displacement  functions.  The  minimization  process  will  minimize  the  gaps 
between  elements  directly  rather  than  through  the  functions  making  them  up. 

The  domain  displacement  field  is  given  by  eqn  (  3.21).  We  use  this  form  exactly  as 
developed  before  and  write  it  as: 

“  1 

-  w  >  =  [I\  [A\  {a}  (3.42) 

.  w  J  n 

where  the  subscript  Q  was  added  to  denote  that  this  is  the  displacement  field  defined  at 
any  point  in  the  domain. 

3.4.2  The  Element  Boundary  Displacements  from  Nodal  Degrees  of 
Freedom 

We  now  establish  an  element  pictured  exactly  as  in  Figure  3.2.  This  time  however,  instead 
of  nodal  degrees  of  freedom  as  defined  in  eqns  (  3.26)-(  3.27)  we  define  them  as 

{^•}T  =  L  uo  v0  <Pi  Vy  w\  i  =  1,2, 3, 4, 5, 6, 7, 8  (3.43) 

and  {A}  is  defined  exactly  as  in  eqn  (  3.28),  with  the  exception  that  now  {A}  is  a  40  term 
vector  rather  than  one  with  44  terms.  This  form  for  the  elemental  degrees  of  freedom  is 
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identical  to  the  standard,  eight  noded  isoparametric  element  used  in  Section  2.3.1.  The 
only  difference  is  that  the  use  of  shape  functions  has  not  been  established.  Instead,  we 
follow  the  method  of  Section  3.3.2  to  define  the  actual  displacements  on  the  boundary 
in  terms  of  the  forty  nodal  degrees  of  freedom.  We  assume  a  quadratic  variation  on  the 
boundary  of  the  element  for  each  of  the  five  degrees  of  freedom.  Now  eqn  (  3.29)  and 
Figure  3.3  can  be  used  to  write 

u0 
Vo 

•fix  '  =  [£MA)  (3-44) 

<Py 

w  J. 

where  again  the  subscript  i  denotes  the  ith  side,  and  the  matrix  [ Q ]  contains  the  necessary 
terms  from  eqn  (  3.29). 

We  cannot  use  eqns  (  3.13)  to  establish  the  displacements  on  the  boundary.  This  is 
because  the  boundary  displacement  functions  in  terms  of  the  nodal  degrees  of  freedom, 
given  in  eqn  (  3.44),  do  not  include  the  w,x  and  w,y  terms.  If  we  were  to  restrict  the 
analysis  to  rectangular  elements  aligned  with  the  global  axis  system,  we  would  have  either 
w,x  or  w,y  on  any  one  side.  This  would  enable  the  displacements  parallel  to  each  side  to 
be  calculated.  A  more  desirable  option  would  be  to  define  new  functions  for  u,  v  and  w. 
Proceeding  in  this  direction,  we  define 

n 

+  53  +  « 

>=i 


(3.45) 
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Wb  =  W0 

where  the  subscript  b  denotes  a  boundary  displacement.  These  expressions  for  the  dis¬ 
placement  fields  are  the  same  as  those  in  eqns  (  3.13)  with  the  w,x  and  w,y  terms  removed, 
making  it  a  C°  continuous  field.  The  displacements  along  the  ith  side  can  now  be  written 
as: 


f  «o  ) 


u  1 

'  1 

0 

ci 

0 

0  ' 

Vo 

V 

»  = 

0 

1 

0 

ci 

0 

< 

<Px 

w  J 

r, 

0 

0 

0 

0 

1 

<Py 

(3.46) 


where  the  subscript  b  has  given  way  to  a  subscript  T,  and  ci  and  ci  are  the  same  as  defined 
eqns  (  3.23)-(  3.24).  Upon  substituting  in  eqn  (  3.44),  the  above  becomes: 


Mr, 


'  1  0  ci  0  0 

0  1  0  ci  0 

0  0  0  0  1 


[S],{A} 


For  convenience,  we  let 


U\ 


1  0  cx  0  0 

0  1  0  ci  0 

0  0  0  0  1 


(3.47) 


(3.48) 


so  eqn  (  3.47)  can  be  written  as: 


3.4.3  The  Least  Squares  Implementation:  Method  II 


(3.49) 


In  the  last  two  sections  we  established  expressions  for  the  displacements:  first,  within 
the  domain  displacement  field ,  (Section  3.4.1),  and  then  for  the  displacement  field  on  the 
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boundary,  (Section  3.4.2).  We  now  proceed  in  the  same  manner  as  in  Section  3.3.3  and 
develop  an  element  domain  displacement  field  in  terms  of  the  nodal  degrees  of  freedom. 
This  time  however,  the  Least  Squares  method  will  be  used  to  minimize  the  difference 
between  the  actual  displacement  fields  on  the  boundary  and  not  the  displacement  functions, 
which  are  their  components. 

We  establish  a  functional,  /,  as  being  the  integral  of  the  squares  of  the  difference  be¬ 
tween  the  domain  displacement  field  evaluated  on  the  element  boundary  and  the  boundary 
di  placement  field  in  terms  of  the  nodal  degrees  of  freedom.  This  can  be  written  as: 

I  =[  [h/2  ({un}r  -  {ur})2  dz  ■  dT  (3.50) 

Jr  J-h/2 


Again,  the  symbol  T  refers  to  boundary  so  the  first  term  in  the  integrand  is  interpreted  as 
the  domain  displacement  field  evaluated  on  the  boundary.  This  equation  is  very  similar  to 
eqn  (  3.34)  with  two  notable  exceptions.  The  first,  as  already  mentioned,  is  the  obvious 
difference  of  the  type  of  variables  in  the  integrand,  actual  displacements  rather  than 
displacement  functions.  The  second  difference,  which  is  a  result  of  the  first,  is  that  the 
integral  now  has  another  dimension  added  to  it.  This  is  because  the  terms  under  the 
integral,  the  displacement  fields ,  are  functions  of  z  as  well  as  x  and  y.  Hence,  in  minimizing 
the  difference  between  the  displacements  on  the  boundary,  the  thickness  of  the  element 
must  also  be  considered.  Upon  substituting  in  eqn  (  3.42)  and  eqn  (  3  49),  we  can  write 
the  above  functional  as: 

1=  f  fk/2  ([l}[A]{a}-[J][G}{A})2dz-dT  (3.51) 

Jr  J-h/2 


Minimizing  this  with  respect  to  the  a’s,  we  get: 
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=  l  J*2  (mt  mT  m  mi  -  iAf  pf  tJi  isi)  •  dr = o  <3.52) 

Again,  the  factor  2  has  been  divided  out.  Also  note  that  the  differential  area  has  been  split 
up  into  a  differential  length  and  thickness  because  it  will  be  convenient  to  perform  the 
thickness  integration  first  to  simplify  the  expression.  Noting  that  the  only  z  dependency 
is  in  [I]  and  \J],  we  can  write: 

J  {[A)T  [M]  [A]  -  [A}t  m  [G])  dr  =  0  (3.53) 

where 

rh/ 2 

[M\  =  /  [X]T  [2]  dz  (3.54) 

J-h/ 2 

Chl2 

[AO  =  /  W?\3\dz  (3.55) 

J-h/2 

Now,  eqn  (  3.53)  can  be  evaluated  in  the  same  manner  as  eqn  (  3.36)  and  the  end  result 
solved  for  the  a's  in  terms  of  the  nodal  degrees  of  freedom.  The  result  is  : 

{a}  =  (M]T[M]M])_1[X]T[A0[a]{A}  (3.56) 

or 

{«}  =  [Hn\  {A}  (3.57) 

where 

\H„  |  =  (Mf  [M]Ml)'1MTlAmC]  (3.58) 


This  result  is  now  the  basis  for  the  second  version  of  the  Least  Square  element.  Just  as 
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with  Method  I,  we  now  have  the  ability  to  write  the  element  displacement  field  anywhere 
in  the  domain  in  terms  the  local  coordinates  and  the  nodal  degrees  of  freedom.  The  result 
is 


(3.50) 


Except  for  the  differences  between  [77/]  and  [H[j\  there  is  no  difference  between  eqn 
(  3.41)  and  eqn  (  3.59).  All  of  the  discussion  at  the  end  of  Section  3.3.3  concerning  the 
implementation  of  the  expressions  for  the  displacements,  is  applicable  to  both  Method  I 
and  II. 

3.5  Preliminary  Numerical  Considerations:  Method  I  -vs- 
Method  II 

The  idea  of  approximating  a  Cl  continuous  element  through  the  Least  Squares  method 
is  the  same  in  both  Method  I  in  Section  3.3  and  Method  II  in  Section  3.4.  However,  the 
details  of  the  developments  are  quite  different,  and  as  would  be  expected,  the  two  methods 
will  result  in  two  different  elements.  It  is  beneficial  to  discuss  some  of  the  differences  at 
this  time. 

First,  and  probably  most  obvious,  the  element  from  Method  I  requires  44  degrees  of 
freedom,  while  the  one  from  Method  II  has  only  40.  Having  fewer  degrees  of  freedom 
may  make  Method  II  computationally  more  efficient.  In  addition,  the  w  displacement  in 
Method  I  is  defined  by  twelve  degrees  of  freedom,  requiring  at  least  twelve  a’s  in  the  initial 
function  for  w.  Method  II  requires  no  more  than  eight  a’s  to  describe  each  unknown, 
as  each  displacement  variable  is  defined  by  eight  degrees  of  freedom.  Requiring  fewer 
unknowns  in  the  initial  representation  of  the  displacement  field  (fewer  a’s)  is  important, 
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for  it  makes  the  matrix  inversion  required  to  find  [Hu]  easier  than  for  that  required  to  find 
[Hj\.  (See  eqn  (  3.37)  and  eqn  (  3.56).)  In  fact,  the  inversion  of  these  matrices  becomes 
a  limiting  factor  into  future  research  into  developing  more  Least  Squares  elements.  In 
eqn  (  3.14)  the  number  of  terms  carried  in  the  expressions  quickly  reaches  a  maximum. 
Since  the  integrations  of  eqn  (  3.36)  and  eqn  (  3.53)  are  carried  out  in  the  local  coordinate 
system,  one  variable  varies  from  -1  to  +1,  while  the  other  is  equal  to  ±1.  Upon  integration, 
the  evaluation  of  odd  terms  go  to  zero,  and  the  even  terms  remain.  Consider  the  matrix 
[O],  where 


[O]  =  [A]t  [A]  (3.60) 

and  [.A]  is  defined  by  eqn  (  3.20).  Upon  integrating  [O]  over  the  boundary  of  the  element, 
as  in  eqn  (  3.36),  it  can  be  shown  that  the  diagonal  is  positive  definite,  because  the  terms 
come  from  the  squares  of  the  individual  terms.  However,  the  off-diagonal  terms  come  from 
the  even  terms  in  the  original  polynomial  expression.  If  n  in  eqn  (  3.14)  is  chosen  so  that 
the  expressions  are  complete  biquadratic3,  then  [O]  is  ill-conditioned,  and  not  invertible. 
The  term  tj2£2  is  the  one  which  causes  this  problem.  So,  in  order  to  get  twelve  unknown 
a’s  in  each  displacement  function  of  eqn  (  3.14),  this  biquadratic  term,  rj2£2,  must  be 
dropped.  The  expression  for  [AJ,  as  shown  in  eqn  (  3.16),  must  be  modified  to  be 

{A}T  =  |AJ  =  U  1  Z  vt  T?  e  r,e  r?3  £3  n\  r,e J  (3.61) 

With  this  form,  [ O ]  is  well  conditioned  can  be  inverted  easily.  This  inversion  problem  is 

by  no  means  new  and  has  previously  been  encountered.  In  an  article  addressing  the  use  of 
5  By  biquadratic  we  mean  the  product  of  complete  quadratic  polynomials  in  each  variable. 
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the  Least  Squares  method  to  smooth  discontinuous  stresses,  Hinton  and  Campbell  (1974 
) [33]  state  that  the  tendency  towards  ill-conditioning  “may  be  overcome  to  some  extent” 
by  using  orthogonal  polynomials  such  as  Legendre  polynomials.  If  formulations  requiring 
a  higher  number  of  alphas  axe  required,  a  few  more  higher  order  odd  terms  may  be 
added  before  [O]  becomes  ill-conditioned.  In  addition,  it  may  be  possible  to  experiment 
with  trigonometric  functions  in  order  to  add  more  unknowns,  if  necessary,  while  still 
maintaining  the  ability  to  invert  the  matrix.  If  judiciously  chosen,  the  trigonometric 
functions  could  be  under  integrated  after  the  inversion  process  to  approximate  the  original 
polynomial  terms. 

Lastly,  but  probably  more  importantly,  we  need  to  consider  the  more  conceptual  dif¬ 
ference  between  the  two  new  Least  Squares  methods.  In  Method  I,  the  Least  Squares  tech¬ 
nique  enforced  the  constituent  functions,  the  displacement  functions  of  the  displacement 
fields,  to  be  compatible  across  element  boundaries.  This  indirectly  enforces  compatibility 
along  element  boundaries  in  the  displacements.  In  contrast,  Method  II  enforces  compati¬ 
bility  directly  to  the  actual  displacement  fields  by  forcing  them  through  the  Least  Squares 
technique  to  fit  a  known  compatible  displacement  field.  Looking  at  the  mechanism  of  the 
methods  more  closely,  we  see  that  Method  I  forces  the  normal  derivative  of  w  along  each 
side  of  the  element  to  be  linear.  In  other  words,  at  the  element  boundary  the  slope  of  the 
plate  will  become  a  linear  function  normal  to  the  edge,  but  it  will  be  quadratic  parallel 
to  it.  Method  II,  on  the  other  hand,  forces  the  displacements  to  be  like  those  with  w,x 
and  w,y=  0  on  the  boundary.  In  effect,  Method  II  will  force  the  w  displacement  to  be 
a  constant  along  the  boundaries  of  the  elements,  if  allowed  to  do  so.  This  sounds  like 
an  unacceptable  approximation.  However,  if  we  choose  w  to  be  an  eight  term  quadratic, 
then  w  cannot  be  constant  on  the  boundaries  except  for  the  trivial  solution  of  w  =  0 
everywhere.  Thus,  through  the  Least.  Squares  method,  the  displacement  fields  represent 
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those  which  minimize  the  gaps  between  the  elements. 

In  light  of  the  above  discussion,  Method  II  is  presented  along  with  Method  I,  not 
for  its  direct  engineering  applicability,  but  as  a  demonstration  of  the  power  of  the  Least 
Squares  method.  Method  II  will  be  shown  to  provide  some  encouraging  results  despite 
the  unacceptable  approximation  discussed  above.  The  intent  is  to  demonstrate  a  second 
method  which  may  be  aoplied  more  effectively  for  other  problems.  It  could  be  greatly 
improved  if  a  better  choice  of  displacement  functions,  based  upon  the  nodal  degrees  of 
freedom,  can  be  found. 

3.6  Finite  Element  Formulation  of  the  Least  Squares  Ele¬ 
ments 

The  development  of  the  stiffness  and  mass  matrices  for  the  Least  Squares  elements  follow 
the  basic  procedure  used  in  Section  2.3.  However,  due  to  the  increase  in  complexity  in  the 
displacement  field,  the  equations  and  matrix  algebra  become  much  more  complicated.  As 
a  result,  the  following  derivation  presents  some  of  the  intermediate  matrices  by  symbol 
only,  and  the  interested  reader  is  referred  to  Appendix  A  for  the  details. 

3.6.1  Stiffness  Matrix  Development 

The  expressions  for  in-plane  strains  within  any  layer  of  an  element  can  be  found  by  sub¬ 
stituting  eqns  (  3.13)  into  the  appropriate  strain  displacement  relations  resulting  in 


k  (  4z3>\ 

4  =  UL  +  J2  ~  ZW'**  +a^r,r  (^2  ~  ^2  J 


-  1l4z 

4  =  Vly  +  Y,p]vly-  zw,yy+PkVy.y  2  ~  ^2 

7=2  \  t 


(3.62) 
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K 

Ixy  =  ul,y  +  vo,x  +  5Z  (PX.y  +  Piv\,x) 


-2zw,xy  +  (a  k<Px,y  +  PkVy 


4'$ 


The  transverse  strains  simply  become 


4-£) 

4-£) 


(3.63) 


In  order  to  keep  future  equations  less  cumbersome,  as  well  as  for  convenience,  the  super¬ 
script  1,  which  implies  the  value  of  the  parameter  for  the  reference  (first)  layer,  will  be 
dropped.  In  other  words,  the  variables  u„,  v0,  <pT  and  <py  will  be  understood  to  refer  to  the 
values  for  the  first  layer  of  the  laminate.  We  can  express  eqn  (  3.62)  in  matrix  notation  as 


(3.64) 


where  the  forms  of  [5]  and  {6}  (see  Appendix  A)  were  purposefully  chosen  to  aid  the 
computations  to  follow.  Similarly,  eqn  (  3.63)  is  written  as 


"  }  =b)‘  =  [Sil  {«'} 


(3.65) 


Referring  to  Appendix  A,  we  see  that  [5]  is  a  3  by  18  element  matrix,  while  [S't]  is  a  2  by 


4  matrix. 
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With  the  above  equations,  the  expression  for  in-plane  strain  energy,  U,  can  be  written 


as: 


V  =  \jvWT{t}dv 

=  \JvUflQ]{‘)dV 

=  \Jv{6}TlS]TmS]{S)dV  (3.66) 


where  [Q]  was  defined  in  Section  2.2.1.  Similarly  the  transverse  strain  energy,  Ut,  becomes: 


V,  =  \fv{y)T[Q]h}dV 

=  \jv  {«'}T  |S|]T  [<3]  [s,|  {«'}  dv  (3.67) 


Next,  we  perform  the  integration  through  the  thickness  of  the  element.  In  doing  so,  the 
results  can  be  written  as: 


U 


\  JA  WT  [«]  w 


dA 


1 

2 


[  {V}7  fn‘l  {V} 

Ja  1  1  1  J  1  J 


dA 


(3.68) 


Ut 


(3.69) 
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where 

rh/2 

[«]  =  /  [S]T[Q][S\dz  (3.70) 

J-h/2 

[fi‘]  =  [o\[St)dz  (3.71) 

The  matrices  defined  in  eqn  (  3.70)  and  eqn  (  3.71)  are  counterparts  to  the  extensional, 
bending  and  coupling  stiffness  matrices  defined  in  eqns  (  2.16)-(  2.19)  in  Section  2.2.1. 
These  [A],  [J3]  and  [D]  matrices  have  specific  physical  interpretations  to  them  and  are 
discussed  in  any  fundamental  composites  textbook  (see  Jones  [41],  Tsai  and  Hahn  [137] 
or  Christensen  [22]).  Unfortunately,  the  components  of  [fl]  and  [fi(]  do  no  lend  themselves 
to  as  nice  a  physical  description.  The  two  matrices  can  be  broken  down  into  extensional, 
bending,  coupling  and  so-forth  submatrices,  but  there  is  nothing  to  be  gained  at  this 
point  from  doing  so.  If  computational  efficiency  for  specific  laminates  were  desired,  then 
knowing  which  submatrices  go  to  zero  for  these  cases  would  be  beneficial.  The  form  of  [5] 
and  [5(]  were  chosen  to  make  the  integrations  in  eqn  (  3.70)  and  eqn  (  3.71)  as  easy  as 
possible. 

The  next  step  in  the  finite  element  formulation  is  to  write  {£}  and  {£f}  in  terms  of 
the  nodal  degrees  of  freedom.  We  first  establish  a  simplified  version  of  {6}  and  {6*}.  We 
define  [Li]  and  [£]]  (see  Appendix  A)  such  that 


{«>  = 

[i.KM 

(3.72) 

M  - 

[l\]  {«:,} 

(3.73) 

{MT  = 


where 
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L  uo,x  uo,y  vo,x  vo,y  Px,z  ^Px,y  Py,z  Py,y  wo,xx  ^o,iy  wo,yy  J  (3.74) 


and 


=  L  Py  U>o,x 


wo,y  J 


(3.75) 


Next,  in  preparation  for  integration  in  the  local  coordinate  system,  the  conversions  of  the 
global  derivatives  to  the  local  coordinate  system  are  established.  We  establish  [£2]  and 
[Z4]  (see  Appendix  A)  such  that 


{^ry}  =  (3.76) 

{4,}  -  [iS]!**}  (3.77) 


where 


{  V)T  = 

L  uo,rj  Uo,£  vo,ri  uo,£  Px,tj  <Pz,£  Py£  Fi(w)  F2(w)  Fs( w)  J  (3.78) 

and 

{^}  =  L  V*  Py  wo, T,  w0,z  J  (3.79) 

In  eqn  (  3.78)  the  terms  Fi(w),  F2{w)  and  Fz(w)  arise  from  the  required  higher  derivatives 
of  to.  They  are  defined  by: 


Fl(w)  —  W,rfri  C3xir]T)  cAVitii ; 
F2(w)  =  w,V£  -c3x,r,{  -cty,^ 
F3(w)  =  w,tf-C3X,tf-C4y,tf 


(3.80) 
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where 

03  =  Tnw,,  +ri2U>,( 

C4  —  ^2l‘Wyrj+T'22 

I\j  =  Componentsof  [J]-1 

These  equations  have  become  complicated  because  of  the  second  derivatives  of  w  in  eqn 
(  3.74).  These  require  new  transformations  relating  the  higher  order  derivatives  in  the 
global  to  the  local  coordinate  systems.  These  transformations  show  up  in  both  eqns 
(  3.80)  and  in  [£2]-  The  complete  details  of  these  transformations,  which  are  not  normally 
seen  in  the  literature  4,  can  be  found  in  Appendix  B.  In  a  final  sequence  into  developing  the 
element  stiffness  matrix,  we  establish  a  new  column  matrix  containing  all  of  the  required 
local  derivatives  of  the  displacement  functions 

=  L  U0,r)  Uo£  Vo,T)  V0j( 

<Px,Z  Py,*)  W1TJ  wyr]T]  J  (3.81) 

This  matrix  is  found  through  differentiation  of  eqn  (  3.40)  with  either  [Hi]  or  [ Hji ]  de¬ 
pending  upon  whether  Method  I  or  Method  II  is  being  used.  The  terminology  [ H]  will 
hereafter  imply  either  Method  I  or  Method  II.  We  can  write 


{«*)  =  [A«]  ["1  I'M 

(3.82) 

K<}  = 

(3.83) 

4 In  the  literature  search  for  this  work,  only  one  article  with  similar  derivations  was  found.  See  Reddy 
(1989)  (116]. 
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where  [Atf]  and  [«4^j  are  differentiated  versions  of  [.4],  and  are  defined  in  Appendix  A. 
The  matrix  {<5^}  is  related  to  {^}  through  the  relation 

{*.*}  =  (3-84) 

where  [£3]  is  given  in  Appendix  A.  Substituting  eqns  (  3.72)-(  3.83)  into  eqn  (  3.68)  and 
eqn  (  3.69),  the  expression  for  the  elements  total  strain  energy,  U3e  =  U  +  Ut,  becomes 

U*'  =  \i  [/*‘f  [Ql]  [?}  dA  (3-85) 

where 

m  =  [LMWVMMW 

M  =  MNKflw 

Minimizing  eqn  (  3.85)  with  respect  to  {A}  gives  the  stiffness  matrix.  The  result  is: 

[*]  =  Ja  ([P]TmP]  +  [pf  [n‘]  [/?*])  dA  (3.86) 


3.6.2  Mass  Matrix  Development 

The  mass  matrix  development  follows  much  the  same  procedure  as  in  Section  2.3.3.  The 
expression  for  kinetic  energy,  Uke,  can  be  written  as 
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Uke  =  ~\u2  f  P  («2  +  *>2  +  u>2)  dV  (3.87) 

or 

Uke  =  ~\U>2  Jv  p{u)T  {“}  dV  (3.88) 

where 

{u}  =  [  u  v  w  J  (3.89) 

Following  the  same  general  procedure  as  in  Section  3.6.1,  we  establish  a  matrix,  [5m],  so 
that  wt  can  write* 

{n}  =  [Sm]{6m}  (3.90) 

where 

{^  }  =  [  U0  V0  W0  1 9X  fly  flX  fly  flx  fly  W,X  W,y  J  (3-9l) 

and  [5m]  is  defined  in  Appendix  A  by  eqn  (  A. 18).  With  these  equations,  eqn  (  3.88) 
becomes: 

uke  =  ~\u2  JvP  {*m}T  [Sm)T  [5m]  {6m}  dV  (3.92) 

This  expression  can  be  integrated  through  the  thickness  resulting  in 

Uke  =  ~^u)2  J^{6m}T  [nm]  {6m}  d A  (3.93) 


where 


rh<2  .  r 

[ftm]  =  /  p[5m]T  [5m]dz 

J-h/2 


(3.94) 
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The  next  task  is  to  represent  {<5m}  in  terms  {A}.  Towards  this  end,  we  establish  a  new 
column  matrix 


$y  W  W„j  W,{  J 


(3.95) 


where 

{^}  =  M][ff]{A}  (3.96) 


In  this  equation,  the  matrix  [.4]  is  exactly  as  defined  in  eqn  (  3.20).  Finally,  {<5m}  is  related 
to  /  <5^  |  through 


(3.97) 


where  [Z/m]  is  defined  in  Appendix  A.  Putting  all  of  the  above  equations  together  allows 
Uke  to  be  expressed  in  terms  of  {A}.  The  result  is 

Uke  =  J JA}t  [/?m]T  [fi m)  [f3m ]  {A}  dA  (3.98) 

where 

\Pm\  =  [Lm]  [A]  [H]  (3.99) 


To  find  the  expression  for  the  mass  matrix,  we  minimize  eqn  (  3.98)  with  respect  to  {A}. 
The  result  is 


M  =  f  [Pm}T  [nm}{Pm}dA 

Ja 


(3.100) 
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This  is  the  final  expression  for  the  element  mass  matrix. 

3.6.3  Stress  Stiffness  Matrix  Development 

The  finite  element  coding  for  Method  I  and  Method  II  will  be  extended  to  include  both 
the  buckling  problem  and  the  pre-stressed  vibration  problem.  To  accomplish  this,  we 
develop  the  initial  stress  stiffness  matrix,  [&*],  following  the  procedure  outline  in  Cook 
[23].  We  assume  that  the  strains  are  composed  of  both  linear  and  nonlinear  portions.  The 
are  written  as: 

{e}  =  Ul]  +  {evi}  (3.101) 


The  nonlinear  portion  of  the  strains  {ejvz;}  includes  the  higher  order  terms  as  in  the 
Lagrangian  strain  definitions.  The  buckling  and  pre-stressed  vibration  analysis  will  be 
limited  to  include  only  the  in-plane  normal  stresses  in  the  following  development,  but 
they  could  easily  be  modified  to  include  in-plane  shear  and  transverse  stresses  quite  easily. 

The  strain  energy  due  to  the  nonlinear  portion  of  eqn  (  3.101)  can  be  expressed  in 
terms  of  the  nonlinear  portion  and  the  pre-stress  present  in  the  element.  Thus 


Unl  = 


(3.102) 


The  nonlinear  portions  of  the  strains  are  written  as 


( «,  1 

L  -  1 

V,x 

w,x 

0 

0 

0 

l  *y  J 

U,  2 

0 

0 

0 

U,y 

V,y 

W,y 
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=  ^[Q]M 


(3.103) 


Prom  eqn  (  3.90)  we  can  write: 


=  [Sm]  ~  } 


(3.104) 


where  {5,™  }  is  the  derivative  of  eqn  (  3.91)  with  respect  to  x.  Similarly  we  can  write: 


which  now  allows  us  write: 


=  [Sml{«,”} 

W,y  I 


(3.105) 


0  U  5 "  1 

J  1  6$  J 


(3.106) 


Next  we  establish 


(3.107) 


where  [/g]  and  }  are  defined  in  Appendix  A.  We  can  now  write: 


1  5m  0 

2  0  5m 


(3.108) 


Next,  in  a  sequence  very  similar  to  that  in  eqns  (  3.73)-(  3.84),  we  establish  the  following 
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relation  to  relate  to  the  nodal  degrees  of  freedom.  (Note:  The  details  of  the  following 

terms  are  omitted,  but  the  matrices  are  very  similar  to  those  in  eqns  (  3.73)-(  3.84)  and 
are  distinguished  through  the  use  of  the  '  symbol.) 

K)  =  [4]  [£'3]  Kf]  W  { A}  (3.109) 

At  this  point  it  is  convenient  to  define: 


W)  =  [In]  {L'2}  [L'3]  [^]  [H]  (3.110) 

so  that  we  can  simply  write 


[/?']  {A}  (3.111) 

. 


This  equation  can  be  substituted  into  eqn  (  3.102)  to  give: 

tto-iXw1 w[t 


We  manipulate  the  last  two  terms  in  the  integrand  to  give: 


["  o xo  0  0  0  0  0 

0  c 7yo  0  0  0  0 

roiT  /  a*°  1  _  0  0  oxo  0  0  0 

1  J  \  Vyo  j  0  0  0  cryo  0  0 

0  0  0  0  gxo  0 

0  0  0  0  0  azo 


[ao]  {lij-y} 


(3.113) 
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Substitution  of  the  expressions  for  {uxy  }  into  eqn  (  3.113)  and  the  result  into  eqn  (  3.112), 
gives  the  required  form  to  provide  a  symmetric  matrix.  The  result  is: 


uNL  =  \{&}T  fvm: 


Sm 

0 

T 

KJ 

Sm 

o 

0 

Sm 

O 

_ 1 

Sm  J 

[?]{*}  W 


(3.114) 


The  element  stress  stiffness  matrix  is  now  established  through  minimization  of  eqn 
(  3.114)  after  integration  through  the  thickness.  The  end  result  is 


[M=  [  W]t[Ob]  W]  dA 

J  A 


(3.115) 


where  the  integration  through  the  thickness  produced  the  matrix  [Ob]  and  is  defined  as 


0  1 
j 


(3.116) 


It  turns  out  that  [Ob]  can  be  written  in  terms  of  the  matrix  [fim]  defined  in  eqn  (  3.94). 
The  relationship  is: 


[Ob]  = 


Pxo  [I^m]  0 

0  Pyo  [flm] 


(3.117) 


Here,  the  terms  pxo  and  py0  are  pre-stress  forces  per  unit  width  stemming  from  the  inte¬ 
gration  of  the  pre-stresses  through  the  thickness.  The  inclusion  of  material  properties  in 
the  displacement  field  has  added  a  complexity  to  the  calculations  involved  in  calculating 
the  stress  stiffness  matrix.  In  most  mechanics  problems  this  matrix  is  found  through  much 
simpler  means. 
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3.7  Finite  Element  Implementation  of  the  Least  Squares 
Elements 

In  the  past  three  sections  we  have  developed  equations  for  the  element  stiffness  (eqn 
(  3.86)),  mass  (eqn  (  3.100))  and  stress  stiffness  (eqn  (  3.115))  matrices.  These  three  ele¬ 
ment  matrices  can  be  assembled  into  the  global  matrices  following  standard  finite  element 
techniques.  The  resulting  equation  to  be  solved  becomes 

([A']  -  w2  [A/]  -  Pa  {A}  =  0  (3.118) 

where  here  we  have  assumed  pzo  =  pyo  =  p0.  If  the  value  of  pa  is  zero,  then  eqn  (  3.118) 
reduces  down  to  the  vibration  eigenvalue  problem  of  eqn  (  2.59)  in  Section  2.4.1.  To 
solve  the  buckling  problem,  we  set  the  second  term  equal  to  zero  and  solve  the  eigenvalue 
problem  for  pa.  For  the  pre-stress  problem,  we  set  p0  equal  to  some  value  les  than  the 
critical  buckling  load  and  solve  the  resulting  eigenvalue  problem  for  w. 

Stress  calculations  for  the  Least  Squares  method  are  done  in  a  manner  similar  to 
that  already  discussed  for  the  Predictor  Corrector  method.  Stresses  are  calculated  at  the 
three  by  three  Gauss  points,  which  are  then  treated  as  the  nodes  for  a  reduced  eight  noded 
element.  The  stress  data  from  these  eight  nodes  are  then  used  to  find  the  derivatives  of  the 
stress.  The  derivatives  of  the  stresses  are  needed  to  put  into  the  equations  of  equilibrium, 
which  are  integrated  to  give  the  transverse  stresses. 

The  formulation  of  the  Least  Squares  element  presented  here  is  unique.  There  exists  no 
formulations  like  it  in  the  literature.  In  fact,  the  only  previous  works  which  were  found, 
used  the  least  squares  process  as  a  tool  to  smooth  finite  element  results,  in  particular 
stresses.  A  few  papers  appeared  over  a  short  period  of  time,  but  none  of  these  implemented 
a  least  squares  approach  to  smoothing  functions  to  approximate  a  continuous  element. 
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Published  works  using  the  least  squares  method  for  smoothing  functions  include  Hinton 
and  Campbell  (1974)  [33],  Lynn  and  Arya  (1973)  [70]  and  (1974)  [69],  Hinton,  Scott  and 
Ricketts  (1975)  [34]  and  Razzaque  (1973)  [109]. 


CHAPTER  IV 


NUMERICAL  RESULTS  AND  METHOD 

VALIDATIONS 

4.1  Study  Approach  and  Preliminary  Information 

The  implementation  of  the  Predictor  Corrector  and  the  two  Least  Squares  techniques  into 
the  finite  element  method  was  presented  in  great  detail  in  Chapters  II  and  III.  This  Chapter 
will  present  the  basic  numerical  findings  of  these  methods.  The  approach  will  be  to 
compare  numerical  results  obtained  with  all  three  of  these  methods  to  existing  data  found 
in  the  current  literature.  The  intent  is  to  provide  an  understanding  into  the  capabilities 
and  limitations  of  the  techniques.  The  primary  study  will  concentrate  on  the  solutions  of 
the  eigenvalue  problems  encountered  when  finding  the  natural  frequency  and  the  critical 
buckling  load.  The  vibration  problem  will  be  analyzed  with  all  three  techniques.  However, 
the  stability  problem  was  only  implemented  into  the  two  Least  Squares  methods,  so  will 
not  be  addressed  by  the  Predictor  Corrector  technique.  In  addition,  there  will  not  be  a 
rigorous  mathematical  convergence  study  but  only  a  demonstration  and  discussion  of  the 
convergence  characteristics. 

In  comparing  the  methods  to  existing  data  found  in  the  literature,  an  attempt  has 
been  made  to  provide  results  for  a  wide  range  of  possible  laminates.  While  the  emphasis 
is  placed  on  thick  laminates,  the  number  of  layers  does  not  necessarily  have  to  be  large. 
The  laminates  studied  range  from  as  little  as  three  layers  to  as  many  as  sixteen.  The 
ply  angles  were  not  limited  to  the  standard  0,  90  and  ±45  degrees,  and  both  symmetric 
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Table  4.1:  Material  Properties. 


MATERIAL 

El/Ej 

Glt/Et 

Gtt/ Et 

VLT 

I'XX 

I 

15 

0.5 

.35 

.3 

.49 

II 

40 

0.6 

.5 

.25 

.25 

III 

25 

0.5 

.2 

.25 

.25 

IV 

4.46 

.566 

.395 

.415 

.49 

V* 

11.49 

.566 

.28 

.38 

.49 

VI 

30 

.6 

.5 

.25 

.25 

VII 

?t 

.6 

.5 

.25 

.25 

♦For  this  material:  Et  =  I.HEt^iv)  and  p  =  0.846p(/v) 
f  Indicates  variable  ratio. 


and  anti-symmetric  laminates  are  considered.  In  addition,  to  a  wide  range  of  material 
properties,  a  laminate  made  up  of  two  separate  materials  (hybrid)  was  also  considered. 
In  choosing  such  a  wide  range  of  variables  for  the  different  laminates,  an  appreciation  for 
the  capability  and  limitations  can  be  realized.  The  material  properties  which  axe  used  are 
listed  in  Table  4.1.  The  ply  layup  of  the  hybrid  laminate  considered  in  this  analysis  is 
[45/-45/0/90/45/-45/0/90]4s.  The  middle  eight  layers  are  Material  IV,  and  the  top  four 
and  bottom  four  layers  are  Material  V.  This  laminate  was  chosen  to  match  that  analyzed 
by  Noor  (1990)  [89]. 

Throughout  this  research  many  different  boundary  conditions  were  considered  and 
investigated.  It  is  worth  mentioning  at  this  point,  with  more  to  follow  later,  that  the 
boundary  conditions  are  believed  to  be  a  source  of  variance  for  the  Least  Squares  methods 
when  comparing  them  to  3 -D  elasticity  solutions.  The  boundary  conditions  used  in  the 
remainder  of  this  work  are  given  in  Appendix  C. 
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4.2  Method  Specific  Behavior  and  Convergence 

4.2.1  Predictor  Corrector  Technique 

The  finite  element  program  for  the  Predictor  Corrector  technique,  as  described  in  Chapter 
II,  was  implemented  in  FORTRAN  on  the  CYBER  990  computer.  Since  the  main  part 
of  the  program  is  based  upon  a  standard  quadratic  isoparametric  Mindlin  plate  element, 
there  is  no  need  to  spend  much  effort  on  evaluating  its  performance.  This  type  of  element 
is  well  understood  and  is  one  of  the  most  widely  used  elements  in  plate  analysis.  Thus, 
we  will  concentrate  on  the  performance  of  the  Predictor  Corrector  method  to  calculate 
accurate  shear  correction  coefficients  and  the  resulting  improvement  in  the  eigenvalues. 

As  discussed  previously  in  Section  2.4.3,  the  shear  correction  coefficients  should  pri¬ 
marily  be  functions  of  the  plate’s  material  properties  and  cross  sectional  shape,  and  hence, 
they  should  be  relatively  constant  throughout  the  plate.  Figure  4.1  provides  some  typical 
examples  of  the  shear  correction  factors  calculated  at  different  locations  in  a  quarter  plate 
model.  The  data  points  represent  a  different  x  or  y  location  in  the  plate.  Data  near 
locations  where  the  shear  stresses  are  zero  are  disregarded,  as  these  are  not  accurate.  One 
can  see  that  the  fines  are  relatively  constant,  and  an  average  value  of  the  shear  correction 
factors  can  easily  be  calculated.  It  was  found  that  it  becomes  difficult  to  calculate  the 
shear  correction  coefficients  anyplace  where  the  transverse  shear  stresses  axe  small.  This 
is  understandable  from  a  numerical  viewpoint  when  one  considers  the  method  involved  in 
calculating  these  numbers.  The  integration  through  the  thickness  of  the  plate  relies  upon 
the  addition  of  many  numbers,  which  results  in  inaccuracies  if  the  numbers  are  small.  The 
coefficients  themselves  are  then  found  by  dividing  a  small  number  by  another.  The  result 
is  likely  to  be  inaccurate.  The  important  point  is  that  the  shear  correction  coefficients 
should  be  calculated  from  the  data  found  surrounding  an  area  in  the  plate  where  the 
stresses  Eire  relatively  high. 
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No  attempt  was  made  to  automate  the  calculation  of  the  shear  correction  coefficients 
into  successive  runs  of  the  finite  element  program.  In  determining  the  natural  frequency 
of  a  particular  laminate,  the  program  was  first  run  with  both  shear  correction  coeffi¬ 
cients  equal  to  one  and  the  resulting  average  shear  correction  coefficients  determined.  The 
program  was  then  rim  a  second  time  with  the  new  coefficients  and  the  updated  natural 
frequency  determined.  The  calculations  involved  in  calculating  the  shear  correction  co¬ 
efficients  were  avoided  in  this  second  run.  It  is  interesting  to  note  that  there  is  no  need 
to  iterate  a  solution  to  find  the  most  accurate  shear  correction  factors.  Noor  (1989)  [88] 
also  found  this  to  be  the  case.  Figure  4.2  presents  the  calculated  transverse  stress  for 
a  laminate  calculated  with  shear  correction  coefficients  equal  to  1.0  and  then  with  the 
coefficients  equal  to  0.79.  Note  that  there  is  no  change  in  the  shape  of  the  curves.  The 
lines  fall  exactly  on  top  of  one  another. 

To  validate  the  program’s  ability  to  calculate  the  shear  correction  factors  in  the  present 
analysis,  the  results  are  compared  to  those  obtained  by  Noor  [88],  For  example,  the  values 
published  by  Noor  for  a  9  layer  crossply  laminate  are  kx  =  0.838  and  ky  —  0.730.  The  data 
from  this  analysis  is  presented  in  Figure  4.1.  The  current  analysis  provides  averages  of 
kx  =  0.841  and  ky  =  0.732.  As  a  further  check,  the  shear  correction  factors  were  calculated 
for  a  homogeneous  aluminum  plate.  These  values  came  out  to  be  very  close  to  5/6. 

The  sensitivity  of  the  natural  frequency  to  the  shear  correction  factor  is  shown  in 
Figure  4.3.  The  figure  shows  that  ujn  is  almost  linear  with  respect  to  the  shear  correction 
factor.  This  is  independent  of  the  number  of  layers.  Based  upon  this  data,  one  could 
update  un  with  relative  confidence  knowing  the  relationship  du >n/dkQ.  The  relationship 
of  the  actual  shear  correction  coefficients  with  the  number  of  layers  is  quite  different, 
however.  Referring  to  Figure  4.4,  one  can  see  that  the  shear  correction  factors  do  not 
follow  any  particular  relationship. 
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It  was  no  surprise  that  the  finite  element  program  converges  rapidly  to  the  first  fun¬ 
damental  frequency.  It  was  found  that  9,  16  and  25  element,  quarter  plate  meshes  yielded 
essentially  the  same  answer.  In  fact,  both  a  nine  element  full  plate  and  quarter  plate 
model  give  the  same  natural  frequency.  The  shear  correction  factors,  however,  do  not 
converge  as  fast,  as  can  be  expected.  Because  of  the  need  for  accurate  data  fitting  in 
calculating  the  coefficients,  at  least  a  16  element  quarter  plate  model  is  needed.  A  finer 
mesh  is  probably  desirable.  A  25  element  mesh  was  used  for  Figure  4.1.  This  fact  presents 
a  problem  for  cases  where  the  ply  layups  preclude  quarter  plate  symmetry,  requiring  then 
full  plate  modeling.  In  order  to  get  accurate  shear  correction  coefficients  at  least  64  or  100 
elements  are  required.  This  was  found  to  be  a  serious  drawback  to  the  Predictor  Corrector 
method.  The  computational  time  and  storage  required  to  solve  a  problem  quickly  became 
large. 

If  this  program  is  to  be  implemented  in  an  analysis,  new  methods  of  curve  fitting 
the  nodal  displacements  should  be  studied  to  provide  more  accuracy  and  flexibility.  The 
present  analysis  was  restricted  to  rectangular  elements  because  of  the  methods  used  in 
curve  fitting  the  data.  There  is  a  lot  of  room  for  improvement  in  calculating  the  shear 
correction  coefficients  and  there  is  no  reason  that  the  method  cannot  be  improved  upon 
significantly.  If  accurate  single  element  shear  correction  coefficient  calculations  can  be 
made,  the  computational  time  can  be  greatly  reduced.  The  need  for  fine  meshes  can  be 
avoided.  For  maximum  efficiency,  multiple  mesh  runs  may  be  required.  A  finer  mesh  may 
be  needed  in  calculating  the  shear  correction  coefficients  and  then  a  coarse  mesh  used  to 
calculate  the  updated  natural  frequency.  Also,  this  technique  could  easily  be  extended  to 
solve  the  general  shell  problem. 
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4.2.2  The  Least  Squares  Method 

The  behaviors  of  the  two  Least  Squares  elements  developed  in  Chapter  III  require  more 
investigation  than  the  standard  isoparametric  Mindlin  element  used  in  the  Predictor  Cor¬ 
rector  program.  The  elements  are  non-conforming  in  the  sense  that  they  axe  not  free 
from  gaps1  between  elements,  but  they  do  differ  in  that  the  gaps  are  minimized  in  a 
least  squares  sense.  In  the  finite  element  method  incompatibilities  between  elements  nor¬ 
mally  result  in  a  reduction  in  stiffness,  but  in  Methods  I  and  II  the  minimization  of  the 
incompatibilities  through  the  Least  Squares  method  should  also  minimize  this  expected 
effect.  However,  this  least  squares  minimization  of  the  gaps  between  the  elements  can  be 
expected  to  change  the  convergence  properties  of  the  elements.  Recalling  the  difference 
between  Methods  I  and  II,  we  can  also  expect  that  the  two  methods  themselves  may  have 
different  convergence  properties.  In  Method  I,  the  normal  derivative  of  the  out-of-plane 
displacement  along  an  edge  of  an  element  was  forced  to  be  linear,  while  the  tangential 
derivative  was  quadratic.  With  a  displacement  field  of  sufficient  order,  these  requirements 
can  be  enforced.  Thus,  the  gaps,  which  are  minimized  by  enforcing  these  requirements, 
should  tend  to  zero.  This  should  be  true  especially  as  the  number  of  elements  increases. 
However,  one  must  also  consider  the  vehicle  used  in  minimizing  the  gaps.  In  minimizing 
the  gaps,  the  element  is  warped  into  a  shape  not  commensurate  with  the  true  solution. 
Thus,  the  element  may  stmt  out  soft  due  to  the  presence  of  gaps,  but  becomes  more  stiff 
as  the  number  of  elements  increases.  It  may  also  converge  to  a  slightly  stiffer  solution. 

Method  II,  on  the  other  hand,  minimizes  the  gaps  by  forcing  the  out-of-plane  displace¬ 
ment  to  be  equal  to  one  where  the  slopes  of  this  displacement  are  zero.  If  the  displacement 
field  is  of  sufficient  order  to  allow  this,  the  slope  of  the  plate  would  tend  to  zero  at  the 
element  boundaries  causing  a  large  increase  in  stiffness.  This  is  indeed  what  happens  and 


'The  word  gap  is  used  loosely  to  refer  to  both  gaps  and  overlaps. 
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is  stated  without  proof.  It  suffices  to  say  that  when  a  high  order  displacement  field  is  used, 
the  element  becomes  unreasonably  stiff  and  gives  very  poor  results.  As  was  discussed  in 
Section  3.5,  if  the  field  for  w  is  restricted  to  being  quadratic,  it  cannot  meet  the  zero  slope 
case  at  all  edges  of  the  element,  and  hence,  through  the  Least  Squares  technique  should 
not  be  too  stiff.  Thus,  Method  II  is  expected  to  have  gaps,  but  the  effect  of  them  should 
be  minimized  through  the  Least  Squares  method. 

The  results  received  from  Methods  I  and  II  follow  the  discussion  in  the  last  two  para¬ 
graphs.  The  gaps  between  elements  are  most  easily  observed  when  considering  thick 
laminates  made  up  of  just  a  few  layers.  Typical  examples  of  the  incompatibilities  for 
Methods  I  and  II  are  given  in  Figure  4.5  and  Figure  4.6.  These  figures,  for  a  three  layer 
cross-ply,  show  that  Method  I  essentially  has  interelement  compatibility,  while  Method  II 
has  definite  gaps  between  the  elements.  One  can  visually  see,  however,  that  the  gaps  are 
indeed  zero  in  a  least  squares  sense.  Figure  4.7  and  Figure  4.12  show  the  convergence  of 
Methods  I  and  II  respectively  for  a  four  layer  cross  ply  laminate.  The  plots  both  show  the 
ratio  of  v2/u%xact  as  a  function  of  E\/ E?-  Note  that  Method  I  produces  a  stiffer  solution 
as  the  mesh  is  refined,  while  Method  II  gets  less  stiff  as  the  number  of  elements  increases. 
In  addition,  note  that  Method  I  is  converging  as  the  mesh  is  refined,  while  Method  II  is 
neither  converging  nor  diverging  to  a  solution  and  is  moving  away  from  the  exact  solution. 
A  different  convergence  study  is  presented  in  Figure  4.9  and  Figure  4.10,  this  time  for  a 
ten  layer  ±45°  laminate.  In  addition,  the  convergence  is  a  function  of  the  ratio  h/b.  In 
Figure  4.9,  Method  I  is  seen  to  converge  as  before  with  mesh  refinement.  In  addition,  it 
is  also  becoming  relatively  stiffer  as  the  thickness  of  the  plate  increases.  We  will  see  later 
that  this  is  a  typical  trend  for  a  simplified  higher  order  theory.  For  thin  plates.  Method  I 
is  shown  to  converge  to  the  exact  solution.  In  Figure  4.10  Method  II  again  exhibits  some 
peculiar  behavior.  This  time  Method  II  is  moving  towards  the  exact  solution,  but  it  shows 
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little  signs  of  convergence.  An  important  point  to  notice  here  is  that  the  performance  of 
the  element  is  beginning  to  look  acceptable  for  the  thinner  plates,  but  it  quickly  diverges 
for  the  higher  ratios  of  h/b.  Finally,  Figure  4.11  and  Figure  4.12  present  yet  another  case. 
This  time  for  the  sixteen  layer  hybrid  laminate.  Note  the  same  trends  are  present  as  in 
Figure  4.9  and  Figure  4.10. 

From  these  observations  we  can  come  to  several  conclusions.  First,  we  can  conclude 
that  Method  I  exhibits  acceptable  convergence  characteristics.  The  formulation  of  the 
element  allows  the  gaps  between  elements  to  be  easily  closed.  However,  because  of  these 
gaps,  the  solution  converges  opposite  to  what  is  normally  found  in  typical  finite  element 
analysis.  Normally  the  finite  element  method  produces  a  solution  guaranteed  to  be  too 
stiff,  but  it  gets  softer  with  mesh  refinement.  The  presence  of  gaps  between  the  Least 
Squares  elements  invalidates  this  guarantee.  Also,  material  anisotropy  has  little  effect  on 
the  accuracy  of  the  solution,  but  the  model  does  appear  to  become  too  stiff  as  the  thickness 
to  width  ratio  of  the  laminate  increases.  Next,  Method  II  only  exhibits  convergence  for  thin 
to  moderately  thick  plates  with  several  layers.  As  the  thickness  to  side  ratio  of  the  plate 
increases,  gaps  between  elements  soften  the  model  significantly.  This  phenomenon  is  also 
compounded  as  the  number  of  layers  decreases.  Despite  this  apparently  poor  capability  of 
Method  II,  further  results  obtained  using  it  will  be  presented  in  the  sections  to  follow.  Its 
performance  in  estimating  natural  frequencies  will  be  seen  to  be  fairly  accurate  for  many 
cases.  It  is  felt  that  Method  II  may  find  use  if  a  better  understanding  of  its  convergence 
can  be  realized.  There  seems  to  be  some  correlation  between  the  mesh  size  and  the  type 
of  problem  being  solved  as  to  maximizing  the  accuracy.  It  is  felt  that  an  acceptable 
element  could  be  realized  if  some  form  of  hp  convergence  is  implemented.  In  other  words, 
if  the  order  of  the  polynomials  used  to  describe  the  displacement  field  is  manipulated 
in  conjunction  with  the  mesh  size,  an  efficient  accurate  element  may  be  developed.  For 
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information  concerning  the  hp- version  of  the  finite  element  method  see  Holzer  et  al  (1990) 
[35]. 

4.3  Comparisons  to  Known  Solutions 

The  literature  covering  the  past  twenty  years  contains  a  lot  of  data  on  the  analysis  of 
thick  laminated  composite  plates.  A  sufficient  amount  of  3 -D  elasticity  data  exists  to  use 
for  comparisons  due  to  the  multitude  of  analysis  techniques  which  have  been  developed 
over  this  period  of  time.  This  3 -D  elasticity  data  is  generally  considered  to  be  the  exact 
solution  to  the  laminate  problem  and  is  available  for  specific  laminates  and  boundary 
conditions.  It  is  with  these  solutions  that  authors  evaluate  the  performance  of  their 
methods.  Unfortunately,  the  majority  of  the  data  are  static  deformation  solutions  for 
various  loading  cases.  This  is  not  to  say  that  there  is  not  a  fair  amount  of  vibration  and 
stability  data  available.  The  only  data  which  is  difficult  to  come  by  are  detailed  stress 
results  based  upon  the  mode  shapes  obtained  from  a  dynamic  analysis.  In  fact,  the  only 
stress  distribution  data  found  is  in  the  form  of  small  plots  for  just  a  few  cases  published 
by  Noor  (1989)  [87,  88],  (1990)  [89]  and  (1973)  [93].  Therefore,  in  the  comparisons  to 
follow,  the  majority  of  the  data  presented  will  be  in  the  form  of  natural  frequencies  and 
critical  buckling  load.  Also,  as  mentioned  above  in  Section  4.1,  the  specific  cases  chosen 
for  comparison  were  picked  to  represent  a  range  of  different  types  of  laminates.  It  by  no 
means  represents  all  that  is  available.  Also,  unless  otherwise  specified,  all  finite  element 
analysis  were  conducted  with  a  quarter  plate  model  with  a  3  x  3  element  mesh. 

We  will  start  by  first  discussing  how  we  Eire  going  to  compare  the  results  of  the  dif¬ 
ferent  techniques  to  exact  as  well  as  to  other  methods.  Figure  4.13  and  Figure  4.14  are 
presented  to  give  a  feel  for  why  direct  comparison  of  the  non-dimensionalized  eigenvalues 
is  not  helpful.  In  these  figures  the  actual  values  for  the  natural  frequency  and  buckling 
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coefficients  axe  so  close  to  the  exact  solutions  that  differences  are  barely  discernible.  To 
overcome  this  problem,  two  different  methods  are  typically  used.  One,  which  was  already 
used  above  in  Section  4.2.2,  is  to  plot  rather  than  the  actual  coefficients.  The 

other  is  to  plot  the  percent  error  relative  to  the  exact  solution.  Both  methods  give  very 
similar  plots,  but  because  the  percent  relative  error  provides  a  well  understood  measure, 
it  will  be  used  in  the  following  comparisons. 

The  percent  relative  error  corresponding  to  ,'igure  4.13  is  presented  in  Figure  4.15. 
This  data  is  for  the  four  layer  symmetric  cross  ply  laminate  with  a  layup  of  the  form 
[0/90/90/0].  The  differences  between  the  different  methods  can  easily  be  seen  in  this 
figure  as  compared  to  the  previous.  In  fact,  the  errors  are  exaggerated.  Included  in  the 
figure  is  data  from  a  simplified  higher  order  theory,  SHOT.  This  method,  developed  by 
Phan  and  Reddy  (1985)  [105],  uses  a  displacement  field  of  the  form  given  in  eqn  (  3.2),  and 
implements  it  into  a  finite  element  program  through  the  use  of  Hermite  cubic  interpolation 
functions.  The  data  included  in  Figure  4.15,  including  the  exact  solution,  comes  from  this 
reference.  From  the  figure  it  can  be  seen  that  Method  II  behaves  quite  poorly,  while 
Method  I,  the  Predictor  Corrector  and  SHOT  methods  all  behave  similarly  with  about  a 
±1%  error.  In  comparing  the  accuracy  of  the  programs  to  predict  the  critical  buckling 
load  in  Figure  4.16,  we  see  something  quite  different.  This  time  both  Method  II  and  the 
SHOT  have  a  3%  to  4%  error,  while  Method  I  is  still  within  a  1%  error. 

The  next  case  considered  is  the  ten  layered  ±45°  laminate.  This  time  we  see  the 
comparison  as  a  function  of  h/a.  The  data  for  this  comparison  was  published  by  Noor 
and  Burton  (1990)  [89].  Included  this  time  are  the  results  for  the  SHOT  as  well  as  those 
for  a  discrete  layer  theory  and  a  simplified  discrete  layer  theory.  As  could  be  expected, 
the  discrete  layer  theory  provides  the  most  accurate  results.  The  simplified  discrete  layer 
theory  performs  rather  poorly,  with  only  one  data  point  on  the  graph.  This  is  because 


101 


the  simplified  discrete  layer  theory,  discussed  by  Noor  [89],  is  assumed  to  be  a  piecewise 
linear  displacement  field  with  no  shear  correction  for  the  individual  layers.  Next,  the  figure 
shows  that  the  Predictor  Corrector  technique,  along  with  Method  II,  remain  within  1% 
of  the  true  solution.  Method  I  and  the  SHOT  approach  a  3%  error  as  the  plate  thickness 
to  width  ratio  becomes  large.  In  all  fairness,  the  improved  accuracy  of  Method  I  over  the 
SHOT  is  not  as  great  in  reality  as  shown  in  the  figure.  This  is  due  to  the  fact  that  the 
mesh  size  used  has  not  allowed  Method  I  to  totally  converge.  (See  Figure  4.9.)  The  same 
can  also  be  said  about  the  surprising  accuracy  of  Method  II.  (See  Figure  4.10.) 

The  last  comparison  is  conducted  on  the  sixteen  layer  hybrid  laminate.  The  results 
are  presented  in  Figure  4.18.  This  figure  is  very  similar  to  the  previous  figure  for  the 
ten  layered  ±45°  laminate.  This  time  however,  all  three  of  the  methods,  the  Predictor 
Corrector,  Method  I  and  Method  II,  all  perform  better  than  the  SHOT  and  the  simplified 
discrete  layer  approach.  These  results  give  credence  to  the  statement  made  in  Section  1.5.1. 
With  the  correct  shear  correction  factor  a  Mindlin  type  deformation  field  can  provide  as 
good,  if  not  better,  results  when  compared  to  a  higher  order  approach.  The  bottom  line 
is  that  both  Method  I  and  the  Predictor  Corrector  approach  can  provide  accuracies  well 
within  a  1%  to  2%  error. 

At  this  point  we  should  consider  known  possible  reasons  why  Methods  I  and  II  may 
vary  from  the  published  exact  solutions.  As  eluded  to  earlier,  the  boundary  conditions 
applied  in  the  Least  Squares  analysis  may  be  a  source  of  some  variances.  Recalling  from 
Section  3.2.2,  the  unknowns  in  the  problem  (the  degrees  of  freedom)  where  chosen  to 
be  referenced  to  the  bottom  layer.  Thus,  when  setting  the  boundary  conditions,  say  for 
instance  setting  u0  =  0,  we  are  fixing  this  variable  equal  to  zero  for  the  bottom  layer  and 
not  for  the  center  of  the  plate,  as  would  normally  be  desired.  Also  setting  <pa  =  0  on 
the  boundary  sets  it  equal  to  zero  in  every  layer  (see  eqns  (  3.9)-(  3.10)).  In  addition, 
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and  probably  most  important,  is  the  fact  that  in  the  process  of  condensing  the  discrete 
layer  displacement  field  down  to  the  simplified  form,  the  variables  u",  u",  <px  and  p™  have 
lost  their  direct  physical  meaning.  For  instance,  U*  is  not  the  membrane  stretching  at 
the  center  of  the  first  layer.  It  is  the  membrane  stretching  of  the  first  layer  evaluated  at 
2  =  0  (the  center  of  the  plate).  Each  of  the  four  variables,  u",  v ",  px  and  p™  ,  within 
a  layer,  represents  a  section  of  a  set  of  curves  defined  thoughout  the  whole  thickness  of 
the  plate.  The  important  point  is  to  realize  that  these  variables  have  different  physical 
meanings  as  a  result  of  the  simplification  process.  In  order  to  prescribe  a  particular 
displacement  boundary  condition,  say  for  u,  one  must  use  the  displacement  fields  and 
solve  the  expression  for  u  to  get  one  parameter,  say  iia ,  in  terms  of  another,  px.  This 
must  be  done  to  get  the  correct  relationship  amongst  the  variables  to  achieve  the  desired 
boundary  condition. 

Another  point,  to  consider  is  the  behavior  of  the  displacement  field  given  in  eqns  (  3.13) 
when  analyzing  antisymmetric  (even  numbered  layers)  laminates.  Looking  at  eqns  (  3.11)- 
(  3.12)  it  can  be  seen  that  for  these  cases,  uQ  and  va  will  have  the  same  value  for  both  layers 
on  either  side  of  the  laminate  centerline.  This  effect  is  not  felt  to  be  significant,  especially 
for  laminates  with  six  layers  or  more.  However,  care  should  be  exercised  when  analyzing 
antisymmetric  laminates  with  just  a  few  layers.  Most  importantly,  stress  analysis  near 
the  center  of  the  plate  should  be  avoided.  This  is  because  at  the  center  of  the  plate  the 
membrane  stresses  may  be  the  predominant  ones. 

4.4  Stress  Calculations 

The  capability  of  an  analysis  method  to  accurately  predict  the  internal  stress  distribution 
through  the  thickness  of  the  laminate  is  almost  as  important  as  its  ability  to  predict  the 
primary  variables,  whatever  they  may  be.  The  intent  of  this  section  is  to  present  some 
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limited  data  as  to  the  ability  of  Method  I  to  provide  accurate  through-the-thickness  stress 
data  for  the  case  of  free  vibration.  Method  II  is  not  considered  in  this  analysis  because  it 
was  found  that  the  eight  term  quadratic  displacement  field  which  was  required  to  achieve 
accuracy  with  this  method  was  not  sufficient  to  provide  accurate  stress  distributions.  It 
was  found  that  Method  II  required  at  least  twelve  terms  in  the  polynomial  expansion 
before  reasonable  stress  distributions  could  be  realized.  However,  as  previously  pointed 
out,  with  this  many  terms  the  model  becomes  too  stiff  and  inaccurate  results  Eire  obtained. 
For  this  reason,  stress  analysis  with  Method  II  will  not  be  discussed  beyond  this  point.  In 
addition,  the  in-plane  stress  analysis  with  the  Predictor  Corrector  technique  is  essentially 
that  of  any  first  order  technique.  Therefore,  aside  from  presenting  some  transverse  stress 
results,  this  method  will  not  be  considered  in  any  detail. 

In  performing  through-the-thickness  stress  e  .s  for  the  free  vibration  problem,  the 
magnitude  of  stresses  can  only  be  lotcrmined  within  a  constant.  Thus,  displacements  and 
stresses  are  generally  prese  ted  after  being  normalized  to  unity  by  dividing  each  variable 
by  its  own  maximum.  This  method  destroys  the  relative  magnitudes  between  the  different 
variables  themselves,  but  in  keeping  with  what  little  data  is  available  in  the  literature  this 
technique  will  be  applied  except  where  noted.  The  through-the-thickness  stress  data  will 
be  presented  for  three  of  the  cases  presented  in  Section  4.3.  The  data  will  be  presented 
in  a  series  of  seven  plots,  namely  u,  v,  an,  <722,  T12,  T13  and  T23,  all  as  a  function  of  the 
z  position  through  the  thickness.  In  all  of  the  following  graphs,  the  horizontal  grid  lines 
correspond  to  the  layer  interfaces. 

Figure  4.19  through  Figure  4.22  present  data  for  the  ten  layer  crossply  laminate.  It 
is  interesting  to  note  that  in  Figure  4.19  the  displacement  field  is  seen  to  be  very  close 
to  linear  for  both  u  and  v.  This  fact  provides  justification  for  using  a  first  order  shear 
theory  to  model  such  laminates.  In  fact,  the  shape  of  the  two  curves  shows  why  the 
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higher  order  and  simplified  higher  order  theories  work  so  well.  Thus,  the  stress  fields 
plotted  in  the  subsequent  figures  vary  little  from  what  would  be  found  using  these  theories. 
The  transverse  shear  plots  in  Figure  4.22  were  calculated  by  integrating  the  equations  of 
equilibrium  as  previously  discussed.  The  transverse  stress  calculations  from  Method  I 
and  the  Predictor  Corrector  technique  are  compared  in  Figure  4.22.  There  are  some 
differences  between  the  two  methods  which  are  felt  to  be  functions  of  the  curve  fitting 
and  differentiation  methods  used  in  the  calculations.  This  particular  case  for  the  ten 
layered  crossply  laminate  happens  to  be  the  only  one  where  vibrational  stress  data  could 
be  accurately  determined  from  graphs.  Data  was  hand  picked  from  enlargements  of  the 
graphs  taken  from  an  article  by  Noor  (1973)  [93].  This  comparison  of  Method  I  to  a  3 -D 
elasticity  solution  is  presented  in  Figure  4.24  for  an  and  in  Figure  4.25  for  t\$.  Note  that 
the  transverse  stress  is  normalized  with  respect  to  the  in-plane  stress  for  this  case.  In 
subsequent  work  by  Noor  [88,  89],  this  is  not  done,  and  data  is  normalized  to  unity.  In 
addition,  the  x-y  coordinates  of  the  location  in  the  plate  where  the  data  is  taken  is  not 
given  in  this  particular  work.  In  this  research  it  was  found  that  this  ratio  was  dependent 
upon  where  in  the  plate  the  data  was  taken.  For  this  reason,  two  horizontal  axis  are 
provided  with  Figure  4.25,  one  corresponding  to  Method  I,  and  the  other  for  the  exact 
solution.  The  scales  have  been  adjusted  to  allow  for  collocation  of  the  maximum  values 
of  T13.  From  these  two  figures,  Method  I  is  seen  to  provide  normalized  stress  values  which 
fall  very  close,  if  not  identically,  on  the  exact  solutions  for  most  locations  throughout  the 
thickness  of  the  plate.  It  is  felt  that  the  method  used  to  calculate  the  stresses  through  the 
thickness  of  the  plate  can  be  improved  upon  through  some  fine  tuning  efforts.  However, 
because  of  the  lack  of  stress  data  to  compare  to,  the  effort  is  not  possible  at  this  time.  The 
force-displacement  finite  element  technique  should  be  applied  for  which  there  is  a  wider 
range  of  data  available. 
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The  next  series  of  graphs  given  in  Figure  4.26  through  Figure  4.30  are  very  similar  to 
the  previous  set  presented  except  this  time  they  are  for  a  symmetric  nine  layer  cross-ply 
laminate.  The  same  observations  and  statements  made  above  apply  to  this  series.  The 
following  set,  however,  is  different.  Figure  4.31  through  Figure  4.35  are  for  the  sixteen 
layer  hybrid  laminate  and  begin  to  show  a  displacement  field  that  is  somewhat  piecewise 
continuous.  This  type  of  displacement  field  can  still  be  modeled  with  a  great  deal  of 
accuracy  with  a  first  order  or  a  higher  order  displacement  field.  Stresses  would  vary 
somewhat  but  not  significantly.  The  real  difference  is  found  in  the  t\ 2  plot  in  Figure  4.33. 
Note  that  the  shear  stress  is  no  longer  a  smooth  function  as  in  the  previous  plots.  This  is 
because  of  the  ±45°  layers  in  this  laminate.  This  phenomenon  emphasizes  the  fact  that  if 
a  laminate  is  to  carry  in-plane  shear,  then  it  is  wise  to  put  in  several  layers  of  other  than 
0°  and  90°  layers.  In  doing  so,  one  will  reduce  large  interlaminar  shear  stresses. 

4.5  Conclusions  From  Method  Evaluation 

The  performance  of  the  Predictor  Corrector  technique  and  the  Least  Squares  methods  is 
very  good.  The  Predictor  Corrector  approach  has  the  capability  to  extend  the  performance 
of  the  simple  Mindlin  isoparametric  element  to  exceed  that  of  higher  order  and  simplified 
higher  order  theories.  The  ability  to  do  so  rests  upon  the  ability  to  calculate  accurate 
shear  correction  coefficients.  In  this  study  a  method  was  developed  to  calculate  accurate 
coefficients  for  the  purpose  of  demonstrating  the  technique  but  was  limited  to  rectangular 
elements.  The  method  also  required  a  finer  mesh  for  accurate  results  than  may  be  required 
using  some  other  techniques.  All  in  all,  the  Predictor  Corrector  technique  has  the  potential 
for  being  an  accurate,  computationally  cheap  method  to  analyze  fiber  reinforced  composite 
laminates. 

The  Least  Squares  methods  also  have  the  ability  to  perform  very  well  in  analyzing 
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composite  laminates.  Method  I  exhibits  better  convergence  characteristics  than  Method 
II,  making  it  the  immediate  choice  to  use.  However,  Method  II  has  the  potential  of  possibly 
being  improved  through  some  form  of  /ip-convergence  finite  element  technique.  If  this  is 
the  case,  its  low  number  of  degrees  of  freedom  give  it  the  potential  to  be  computationally 
attractive.  Method  I  of  the  Least  Squares  technique  can  give  very  accurate  results,  very 
efficiently,  for  composite  laminates.  Its  performance  is  slightly  better  than  the  higher 
order  or  simplified  higher  order  approaches. 

Of  particular  significance  is  the  ability  of  the  Least  Squares  method  to  approximate 
a  C1  continuous  element.  The  Least  Squares  element  has  many  advantages  over  current 
methods  used  to  solve  the  C1  continuity  problem,  and  it  can  be  extended  to  other  ap¬ 
plicable  areas  where  this  problem  exists.  The  Least  Squares  element  is  not  in  any  way 
limited  to  the  composite  laminate  problem. 


CHAPTER  V 


NATURAL  FREQUENCY  AND  STABILITY 
STUDIES  FOR  LAMINATED  COMPOSITE 

PLATES 

5.1  Study  Approach 

The  effect  of  optimizing  the  behavior  of  a  laminated  fiber  reinforced  composite  plate  is  an 
issue  which  has  not  been  addressed  in  detail  over  the  past  twenty  years.  The  majority  of 
the  work  in  this  area,  relative  to  vibration  and  stability,  has  been  to  show  simply  what 
effect  ply  angle  has  on  natural  frequency  or  buckling  load  for  different  a/b  ratios.  It  is  well 
accepted  that  for  a  square  simply  supported  plate,  the  natural  frequency  is  maximized  for 
9  =  ±45°.  If  the  aspect  ratio  of  the  plate  changes,  so  does  the  value  of  the  angle  which 
maximizes  the  frequency.  This  behavior  is  accepted  and  makes  logical  sense.  Studies 
of  this  type  have  been  published  by  Whitney  and  Leissa  (1969)  [143],  Bert  (1977)  [12], 
Khdier  (1988)  [50]  and  Grenestedt  (1989)  [29]  to  name  a  few.  It  should  be  stated  that 
some  of  these  studies  considered  only  thin  plate  theory.  This  phenomenon  will  be  seen 
in  some  of  the  data  to  follow.  In  addition,  a  typical  study  found  in  the  literature  will 
consider  the  frequency  not  only  as  a  function  of  9  and  a/b,  but  it  will  also  vary  the  plate 
thickness  ratio,  a/h ,  and  ratio  of  anisotropy,  E\/Ei-  Also,  Liew  et  ai  (9189)  [64]  have 
presented  similar  studies  for  triangular  plates.  Some  studies  have  limited  the  number  of 
variables  which  were  considered.  For  instance,  Jones  (1973)  [40]  studied  only  the  effect  of 
a/b  and  the  number  of  layers  on  vibration  and  buckling  of  cross  ply  laminates.  Buckling 
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coefficient  studies  are  not  as  common  as  vibration  studies  but  are  available.  Reddy  and 
Phan  (1985)  [117]  and  Whitney  and  Leissa  (1969  [143]  present  graphs  of  critical  buckling 
load  as  functions  of  0  and  a/b,  while  Noor  (1974)  [94]  presents  buckling  coefficients  as 
functions  of  E\/E2,  h/b,  as  well  as  for  various  values  of  a/b.  Studies  including  pre-stress 
effects  on  plate  natural  frequency  are  the  least  common.  Dawe  and  Craig  (1986)  [25] 
looked  at  pre-stress  effects  on  natural  frequency  while  varying  thickness  and  aspect  ratios. 
Chelladurai  et  al  (1984)  [20]  have  looked  at  fiber  orientation  and  pre-stress  for  different 
aspect  ratios,  but  did  so  only  for  single  layer  lamina.  Studies  of  these  types  are  by  no 
means  limited  to  the  ones  mentioned  above.  Many  more  exist,  and  the  ones  listed  are 
intended  to  provide  some  basic  references  to  the  work  that  has  been  done. 

One  can  see  from  the  above  discussion  that  to  study  the  optimization  of  the  funda¬ 
mental  frequency  and/or  buckling  load  is  not  a  simple  prospect.  There  are  a  wide  range 
of  variables  involved.  One  could  consider  ply  angle,  stacking  sequence,  material  prop¬ 
erties,  thickness  ratios,  aspect  ratios,  boundary  conditions  and  different  combinations  of 
pre-stress.  A  full  mathematical  optimization  study  would  be  insurmountable.  At  best 
we  can  only  try  to  understand  how  these  variables  each  effects  the  behavior  of  the  plate. 
The  design  engineer,  with  a  basic  understanding  of  how  they  all  influence  laminate  be¬ 
havior,  can  then  begin  to  optimize  his  design.  The  intent  of  this  chapter  is  to  present 
some  parametric  studies  to  establish  some  new  understanding  into  how  some  of  these  pa¬ 
rameters  effect  the  vibration  and  stability  of  laminated  composite  plates.  The  study  will 
not  attempt  to  reproduce  findings  already  available  in  current  literature,  but  it  will  try  to 
provide  some  new  insight  into  laminate  behavior. 

For  the  following  study  we  will  hold  the  material  properties,  thickness  ratio  and  number 
of  layers  as  constants.  In  other  words,  we  will  study  a  laminate  consisting  of  a  fixed  number 
of  layers  with  the  material  properties  and  plate  thickness  ratio  also  being  held  as  constant. 
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Table  5.1:  Ply  Angle  Analysis  Data/Figure*  Correlation  Grid. 


CASE 

BNDRY 

COND 

FREE 

VIB 

PRE¬ 

STRESSED 

VIBt 

UNIAXIAL 
BUCKLING 
Ny  =  0 

BIAXIAL 
BUCKLING 
Nx  =  Ny 

A 

a/b  =  0.7 

6 

5.1 

5.2 

5.1 

5.1 

13 

5.3 

5.4 

5.3 

5.3 

B 

a/b  =  1.0 

6 

5.5 

5.6 

5.5 

5.5 

13 

5.7 

5.8 

5.7 

5.7 

C 

a/b  =  1.43 

6 

5.9 

5.10 

5.9 

5.9 

13 

5.11 

5.12 

5.11 

5.11 

D 

a/b  =  1.7 

6 

5.13 

5.14 

5.13 

5.13 

13 

5.15 

5.16 

5.15 

5.15 

♦Numbers  indicate  Figure  number  data  is  displayed  in. 

fSee  Table  5.2  for  the  definitions  of  the  four  pre-stressed  vibration  cases. 


The  aspect  ratio  will  vary  along  with  the  ply  angle  and  boundary  conditions.  Specifically, 
the  study  will  look  at  a  six  layered  laminate  with  the  properties  of  Material  II.  The 
thickness  ratio,  b/h,  was  chosen  to  be  10.  The  layup  is  established  as  [+9/  —  9/  +  9 . . .], 
where  9  will  be  varied.  The  aspect  ratios  considered  will  be  a/b  =  0.7,  1.0,  1.4286  and 
1.7.  Boundary  conditions  considered  were  simply  supported  and  clamped.  (All  four  sides 
the  same.)  Data  was  calculated  for  the  free  vibration  natural  frequencies,  critical  uniaxial 
and  biaxial  buckling  loads  and  natural  frequencies  for  four  pre-stress  cases.  The  pre-stress 
cases  were  defined  relative  to  the  lowest  uniaxial  buckling  load  for  each  case.  Two  pre¬ 
stress  loads,  one  50%  and  one  85%  of  the  lowest  uniaxial  buckling  load,  were  used  for 
both  the  uniaxial  and  biaxial  pre-stress  conditions  giving  the  four  pre-stress  cases  (see 
Table  5.2).  A  tabular  form  of  the  data  collected  is  presented  in  Table  5.1.  The  table  also 
provides  the  number  of  the  figure  on  which  the  data  appears. 
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Table  5.2:  Ply  Angle  Analysis  Pre-stress  Conditions. 


NOMENCLATURE 

Nx 

Ny 

PS1A 

0.50  Nmin 

0 

PS1B 

0.85  Nmin 

0 

PS2A 

0.50  Nmin 

0.50  Nmin 

PS2B 

0.85  Nmin 

0.85  Nmin 

Note:  Nmin  indicates  the  lowest  uniaxial  buckling  load. 


5.2  Initial  Data  Trends 


All  of  the  results  for  cases  A,  B,  C  and  D  are  presented  in  the  next  series  of  sixteen  figures. 
Four  figures  are  provided  for  each  aspect  ratio.  From  Figure  5.1  and  Figure  5.2  we  see  that 
the  natural  frequency  for  free  vibration  and  all  of  the  pre-stress  cases,  as  well  as  the  uniaxial 
and  biaxial  buckling  loads,  are  all  maximized  at  a  value  of  0  equal  to  about  32°.  This 
is  for  Case  A  and  the  simply  supported  boundary  condition.  For  the  clamped  boundary 
condition,  Figure  5.3  and  Figure  5.4  show  that,  even  though  the  natural  frequency  is 
maximized  at  this  same  point,  the  two  buckling  cases  have  shifted.  The  uniaxial  buckling 
load  is  maximized  at  about  28°,  while  the  biaxial  buckling  load  is  maximized  around  39°. 
The  data  for  the  square  laminate,  case  B,  is  presented  in  Figure  5.5  through  Figure  5.8. 
Note  similar  trends  here,  except  that  the  natural  frequencies  and  biaxial  buckling  load  will 
be  maximized  at  45°.  One  other  interesting  observation  will  be  made  at  this  time  and  will 
continue  to  be  observed  in  future  graphs.  We  see  that  in  Figure  5.6  the  natural  frequency 
curves  are  bell  shaped  curves  for  the  50%  buckling  loads  but  are  parabolic  for  the  85% 
cases.  This  phenomenon  even  causes  two  of  the  curves  to  cross  over  one  another  for  small 
and  large  values  of  6.  This  will  not  be  discussed  here  but  merely  noted.  Moving  on  to 
Figure  5.9  and  Figure  5.10,  we  now  see  that  the  point  at  which  the  uniaxial  buckling  load 
is  maximized  is  still  around  30°,  while  the  biaxial  buckling  load  and  free  vibration  natural 
frequency  are  maximized  between  50°  and  60°.  The  pre-stressed  vibration  cases  are  all 
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maximized  around  60°.  For  the  clamped  boundary  conditions  Figure  5.11  shows  similar 
results  with  the  exception  that  the  natural  frequency  does  not  fall  off  at  the  higher  ply 
angles  but  continues  to  rise.  This  increase  is  negligible,  however.  Figure  5.12  shows  that 
the  pre-stress  cases  are  all  maximized  around  67°.  Again,  notice  in  this  figure  the  sharp 
drop  off  in  the  natural  frequency  for  the  two  high  pre-stress  cases.  Case  D,  presented  in 
Figure  5.13  through  Figure  5.16,  shows  the  same  trends  as  found  for  Case  C. 

5.3  Effects  of  Plate  Aspect  Ratio 

We  can  plot  all  of  the  data  presented  in  Figure  5.1  through  Figure  5.16  in  such  a  way 
as  to  more  clearly  show  the  effect  of  plate  aspect  ratio.  For  instance,  Figure  5.17  shows 
the  free  vibration  frequencies  for  the  simply  supported  cases.  From  this  figure  we  can 
clearly  see  how,  as  the  aspect  ratio  of  the  plate  increases,  the  optimum  9  moves  towards 
the  right.  This  is  the  trend  discussed  earlier.  This  trend  is  also  visible  in  Figure  5.18 
for  the  clamped  cases.  Figure  5.19  shows  an  interesting  result.  The  optimum  9  for  the 
uniaxial  buckling  load  does  not  change  with  aspect  ratio  for  the  simply  supported  cases. 
In  addition,  after  the  maximum  buckling  load  is  achieved,  the  buckling  load  is  the  same 
for  all  aspect  ratios.  For  the  clamped  plates,  Figure  5.20  shows  that  the  optimum  theta 
does  vary  slightly,  as  do  the  buckling  loads  past  the  optimum  9.  For  the  biaxial  buckling 
loads,  Figure  5.21  shows  that  the  optimum  value  of  9  follows  trends  similar  to  the  natural 
frequency,  as  shown  in  Figure  5.17.  It  is  interesting  to  note  that  after  the  optimum  theta 
is  achieved,  the  buckling  loads  converge  to  the  same  values  regardless  of  the  aspect  ratio. 

5.4  Effects  of  Pre-Stress  on  Fundamental  Frequency 

The  effects  of  pre-stress  on  the  fundamental  frequency  of  simply  supported  and  clamped 
plates  can  be  seen  in  the  previous  figures.  Looking  back,  Figure  5.8  represents  a  typical 


112 


example  of  such  a  plot.  It  is  no  surprise  that  the  natural  frequency  of  the  plate  decreases 
with  increased  pre-stress  load.  It  is  of  interest,  however,  that  the  shape  of  the  curve 
changes  as  the  pre-stress  load  is  increased.  As  mentioned  previously,  for  the  larger  pre¬ 
stress  loads  the  curve  no  longer  maintains  a  bell  shaped  curve  but  more  parabolic,  and  it 
drops  off  rapidly  at  the  ends.  This  trend  becomes  important  when  considering  the  case 
given  in  Figure  5.16.  If  a  design  engineer  wants  to  choose  the  optimum  9  to  maximize 
the  free  vibration  natural  frequency  for  this  particular  aspect  ratio  plate,  he  would  be  led 
to  pick  9  =  90°.  If  the  plate  becomes  significantly  loaded  either  uniaxially  or  biaxially, 
then  serious  degradation  in  the  vibration  characteristics  would  occur.  Clearly  a  better 
choice  of  9  would  be  65°  for  this  case.  The  important  point  is  the  following:  In  optimizing 
natural  frequency,  one  should  consider  the  loading  conditions  which  will  be  present  on  the 
laminate  in  operation,  as  well  as  the  boundary  conditions. 

5.5  Optimization  of  Fundamental  Frequency  and  Buckling 
Loads 

From  the  above  study  one  can  see  that  optimizing  the  vibrational  and  stability  charac¬ 
teristics  of  a  laminated  composite  can  become  a  complicated  task.  In  fact,  optimization 
of  one  parameter  may  lead  to  poor  performance  in  another.  For  instance,  consider  the 
simply  supported  condition  for  Case  D,  presented  in  Figure  5.13.  If  the  ply  lay-up  angle 
is  chosen  to  be  60°,  then  both  natural  frequency  and  biaxial  buckling  loads  are  very  near 
their  maximum.  However,  the  uniaxial  buckling  capability  of  the  plate  has  been  reduced 
by  approximately  43%.  If  we  choose  a  ply  lay-up  angle  of  30°,  then  the  uniaxial  buckling 
capability  is  maximized,  while  the  fundamental  frequency  and  biaxial  buckling  load  have 
been  degraded  by  21%  and  38%  respectively. 

From  the  above  analysis  we  conclude  that  the  60°  lay-up  angle  is  best  from  a  vibration 
and  biaxial  buckling  point,  while  a  30°  angle  maximizes  the  uniaxial  buckling  capability. 
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One  may  wonder  if  better  overall  performance  could  be  achieved  if  both  60°  and  30°  plies 
were  used  in  the  laminate.  If  we  consider  a  laminate  of  the  form  [—60/  +  30/  —  30)45,  we 
find  the  non-dimensional  natural  frequency  becomes  12.6,  while  the  uniaxial  and  biaxial 
coefficients  become  27.1  and  12.1  respectively.  These  numbers  translate  to  a  10%  reduction 
in  natural  frequency,  a  19%  reduction  in  biaxial  buckling  and  a  19%  reduction  in  uniaxial 
buckling  capability.  Another  option  would  be  to  consider  a  laminate  of  the  form  [—30/  + 
60/  —  60)45.  For  this  lay-up  we  find  the  natural  frequency  and  the  uniaxial  and  biaxial 
coefficients  are  11.8,  27.7  and  10.7  respectively,  corresponding  to  15.6%,  17.3%  and  11.8% 
reductions  from  their  individual  maximums.  For  a  final  case,  consider  a  twelve  layer 
laminate  of  the  same  thickness.  The  lay-up  is  chosen  to  be  [—60/  —  30/  +  60/  +  30/  — 
60/  —  30)45.  This  tim'-  1  i  find  a  non-dimensional  natural  frequency  of  12.9  and  buckling 
coefficients  equ' .  *•  29.3  for  the  biaxial  case  and  11.6  for  the  uniaxial.  These  numbers 
correspond  to  7.7%,  12.4%  and  21.4%  reductions  in  capabilities.  This  last  case,  however, 
has  introduced  a  new  variable,  the  number  of  layers,  which  we  have  purposefully  avoided. 

The  specific  numbers  and  percentages  in  the  above  crude  analysis  are  not  meant  to 
provide  hard  and  fast  numbers,  but  they  are  intended  to  provide  insight  into  how  the 
natural  frequency  and  critical  buckling  loads  can  be  affected  by  ply  angle  and  stacking 
sequence.  The  design  engineer  has  countless  combinations  of  these  and  other  parameters  to 
consider.  In  designing  a  pressure  vessel,  a  biaxial  state  stress  may  be  such  that  Nx  =  2 Ny. 
For  such  a  case  the  optimum  design  would  be  different  yet.  One  important  point  which 
must  be  made  here  is  that  the  behaviors  observed  in  Figure  5.1  through  Figure  5.22  are  for 
a  laminate  with  a  specific  number  of  layers,  thickness  ratio  and  set  of  material  properties. 
If  any  of  these  are  changed,  the  trends  established  above  could  either  be  eliminated  or 
accentuated.  In  addition  new  trends  could  be  observed. 
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5.6  General  Observations  and  Conclusions 

After  investigating  the  behavior  of  laminated  composite  plates  for  many  different  cases, 
some  not  presented  in  this  research,  several  conclusions  have  been  reached  concerning  op¬ 
timizing  design.  First  of  all,  the  optimization  of  the  performance  of  a  laminated  composite 
plate  is  not  a  simple  process.  Some  basic  rules  always  apply.  To  maximize  the  natural 
frequency  or  buckling  load  of  specific  laminate  for  a  given  aspect  ratio  and  set  of  bound¬ 
ary  conditions,  some  ply  angle  will  be  optimum.  The  more  plies  which  are  at  this  angle, 
the  better.  If  there  are  a  fixed  number,  and  the  laminate  contains  other  angles  also,  the 
frequency  is  increased  if  the  optimum  angle  plies  are  moved  towards  the  top  and  bottom 
of  the  laminate.  Next,  the  optimum  angle  for  the  lamina  is  dependant  upon  many  factors 
including  boundary  conditions,  aspect  ratio  and  number  of  layers.  To  optimize  a  specific 
laminate,  the  specific  conditions  under  which  it  will  be  subjected  must  be  considered.  Also 
of  great  importance,  the  optimum  ply  angles  and  stacking  sequence  for  maximizing  the 
natural  frequency  will  not  be  the  same  required  to  maximize  the  buckling  loads.  Finally, 
pre-stress  can  have  a  significant  effect  on  the  frequency  of  a  laminated  plate,  and  the  effect 
must  be  investigated  for  the  specific  case. 


CHAPTER  VI 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  work  conducted  for  this  research  has  shed  new  light  upon  the  effectiveness  of  different 
analysis  methods  used  on  fiber  reinforced  composite  plates.  The  trend  towards  increased 
accuracy  is  driving  the  analysis  methods  to  more  computationally  intensive  approaches. 
This  need  not  be  the  case,  especially  in  the  area  of  thick  composite  laminates.  The 
Predictor  Corrector  technique,  implemented  together  with  the  finite  element  method  and 
the  Least  Squares  elements,  are  just  two  ways  in  which  accurate  results  can  be  realized  for 
little  more  effort  than  for  a  simple  Mindlin  plate  element.  Results  obtained  using  these 
methods  can  be  every  bit  as  accurate  as  techniques  with  increased  complexity  and  many 
more  degrees  of  freedom. 

Of  major  importance,  the  Least  Squares  elements  are  not  limited  to  analyzing  lam¬ 
inated  composite  plates.  The  Least  Squares  elements  developed  in  this  work  have  been 
shown  to  be  a  numerically  effective  method  to  approximate  a  C1  continuous  element. 
This  is  an  important  contribution  to  the  finite  element  field  and  can  be  applied  in  any 
situation  where  C1  continuity  is  required.  In  reality,  the  Least  Squares  technique  could  be 
extended  to  even  higher  orders  of  continuity  with  the  proper  choice  of  functions  describing 
the  primary  variables. 

The  work  performed  here  in  developing  the  Least  Squares  methods  needs  to  be  ex¬ 
tended  in  further  research.  The  element  from  Method  I  should  be  immediately  imple¬ 
mented  into  a  static  force-displacement  finite  element  program.  In  doing  so,  one  would 
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also  develop  moment  curvature  relations  and  plate  constituative  equations.  Such  relation¬ 
ships  would  be  necessary  to  implement  force  and  moment  boundary  conditions,  and  to 
interpret  force  results.  With  the  increased  data  available  in  the  literature,  better  through- 
the-thickness  stress  comparisons  can  be  made  allowing  fine  tuning  to  be  done  to  the  stress 
calculations.  Doing  this  would  also  provide  more  validation  data  for  the  technique.  Mod¬ 
ification  of  Method  II,  through  the  use  of  an  hp-convergence  technique,  should  also  be 
investigated.  Finally,  a  more  thorough  convergence  study  should  be  conducted  to  fully 
understand  the  Least  Squares  method. 

Also  included  in  this  work  was  data  showing  the  relationships  of  ply  lay-up  angle, 
boundary  conditions  and  plate  aspect  ratio  on  natural  frequency  and  buckling  loads.  The 
effect  of  pre-stress  on  the  natural  frequency  was  also  included.  Several  interesting  behav¬ 
iors  were  documented  but  are  restricted  to  a  specific  thickness  ratio,  number  of  layers 
and  type  of  material.  More  work  of  this  type  needs  to  be  conducted  especially  in  areas 
not  common  in  the  literature.  More  studies  into  the  effects  of  boundary  conditions,  pre¬ 
stress  and  simultaneous  optimization  of  natural  frequency  and  buckling  loads  must  be 
conducted.  In  addition,  all  of  these  studies  should  begin  to  look  at  the  stress  distributions 
through  the  thickness  of  the  laminates  so  that  propensities  towards  delamination  can  be 
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Figure  4.1:  Shear  Correction  Coefficient  Variation  Across  Plate.  Material  I  a/b 
layer,  [90/0/90/0/ . . .],  h/b  =  0.2,  BC-1,  5x5  quarter  plate  model. 
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Figure  4.2:  Comparison  of  Transverse  Shear  Stress  with  Different  Shear  Correction  Fac¬ 
tors.  Material  I  a/6  =  1,  10  layer,  [90/0/90/0/...],  h/b  =  0.2,  BC-1,  3x3  quarter  plate 
model. 
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Figure  4.3:  Sensitivity  of  Natural  Frequency  to  Shear  Correction  Coefficient.  Cross  ply 
laminate,  Material  II  a/6  =  1,  h/b  =  0.2,  BC-1,  3x3  quarter  plate  model. 
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Figure  4.5:  Method  I  Interelement  Compatibility.  [0/90/0]  simply  supported  a/6  = 
Material  III,  b/h=4,  Location:  x  =  0.1G6G7,  y  =  0.133333,  3  x  3  quarter  plate  model. 
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Figure  4.6:  Method  II  Interelement  Compatibility.  [0/90/0]  simply  supported  a/b  = 
Material  III,  b/h=4,  Location:  x  =  0.1C667,  y  =  0.133333,  3  x  3  quarter  plate  model. 


Figure  4.8:  Method  II  Convergence.  [0/90/90/0]  simply  supported  (BC-1),  a/b 
Material  VII,  b/h=5,  quarter  plate  model. 


125 


h/b 


Figure  4.9:  Method  I  Convergence.  Ten  layer  [-45/  4-  45/  —  45/...]  simply  supported 
(BC-1),  a/b  =  1,  Material  I,  quarter  plate  model. 
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Figure  4.10:  Method  II  Convergence.  Ten  layer  [-45/  +  45/  -  45/ . . .]  simply  supported 
(BC-1),  a/b  =  1,  Material  I,  quarter  plate  model. 
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Figure  4.14:  Method  Comparison  of  Non-dimensional  Buckling  Coefficient  -vs-  Anisotropy 
Ratio.  [0/90/90/0],  Material  VII,  b/h  =  5,  Aj  =  n\b2/h3E2,  a/b  =  1,  Simply  supported 
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Figure  4.17:  Method  Comparison  of  Percent  Error  in  Natural  Frequency  -vs-  Thickness 
Ratio.  10  layer,  [—45/  +  45/  —  45/ . . .],  Material  I,  a/b  =  1,  Simply  supported  (BC-1) 


Figure  4.18:  Method  Comparison  of  Percent  Error  in  Natural  Frequency  -vs-  Thickness 
Ratio.  16  layer  Hybrid-  Materials  IV  &  V,  a/b  =  1,  Simply  supported  (BC-1) 
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Figure  4.19:  In-plane  Displacements  -vs-  z  Location.  10  layer,  Material  I,  [90/0/90/. . .], 
a/b  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  i  =  0.29167, 


y  =  0.29167. 
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Figure  4.20:  In-plane  Stresses  -vs-  z  Location.  10  layer.  Material  I,  [90/0/90/  . . .],  a/b  =  1, 
h/b  =  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  =  0.29167,  y  =  0.29167. 
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Figure  4.21:  In-plane  Shear  stress  -vs-  z  Location.  10  layer,  Material  I,  [90/0/90/...], 
a/b  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  =  0.29167, 
y  =  0.29167. 
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Figure  4.22:  Transverse  Shear  Stresses  -vs-  z  Location.  10  layer,  Material  I,  [90/0/90/  . . .], 


a/b  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  =  0.29167, 


y  =  0.29167. 


139 


t  /|t  | 

13/  1  1 3lMax 


t  /It  I 

237  1  23  'Max 


Figure  4.23:  Method  Comparison  of  Transverse  Shear  Stresses.  10  layer,  Material  I, 
[90/0/90/  . . .],  a/6  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model, 
1  =  0.29167,  y  =  0.29167. 
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Figure  4.24:  In-piane  Stress  Comparison:  Method  I  -vs-  Elasticity.  10  layer,  Material  VI, 
[90/0/90/ . . .],  a/b  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  4x4  quarter  plate  model, 
x  =  0.29167,  y  =  0.29167. 
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Figure  4.25:  Transverse  Shear  Stress  Comparison:  Method  I  -vs-  Elasticity.  10  layer, 
Material  VT,  [90/0/90/ . . .],  a/6  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  4x4  quarter 
plate  model,  z  =  0.29167,  y  =  0.29167. 
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Figure  4.26:  In-plane  Displacements  -vs-  z  Location.  9  layer,  Material  I,  [90/0/90/ . . .], 
a/6  =  1,  h/b  —  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  =  0.29167, 
y  =  0.29167. 
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Figure  4.28:  In-plane  Shear  stress  -vs-  z  Location.  9  layer,  Material  I,  [90/0/90/...], 
a/b  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  —  0.29167, 
y  =  0.29167. 
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Figure  4.29:  Transverse  Shear  Stresses  -vs-  z  Location.  9  layer,  Material  I,  [90/0/90/  . . .], 
a/6  =  1,  h/b  —  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  =  0.29167, 
y  =  0.29167. 
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Figure  4.30:  Method  Comparison  of  Transverse  Shear  Stresses.  9  layer,  Material  I, 
[90/0/90/  ...],  a/6  =  1,  h/b  =  0.2,  Simply  supported  (BC-1),  3x3  quarter  plate  model, 
x  =  0.29167,  y  =  0.29167. 
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Figure  4.32:  In-plane  Stresses  -vs-  z  Location.  16  Layer  Hybrid,  Materials  IV  V, 
a/6  =  1,  h/b  =  0.3,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  =  0.29167, 
y  =  0.29167. 
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Figure  4.33:  In-plane  Shear  stress  -vs-  z  Location.  16  Layer  Hybrid,  Materials  IV  &  V, 
a/6  =  1,  h/b  =  0.3,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  i  =  0.29167, 
y  =  0.29167. 
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Figure  4.34:  Transverse  Shear  Stresses  -vs-  z  Location.  16  Layer  Hybrid,  Materials  IV  & 
V,  a/6  =  1,  h/b  =  0.3,  Simply  supported  (BC-1),  3x3  quarter  plate  model,  x  =  0.29167, 
y  =  0.29167. 
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Figure  4.35:  Method  Comparison  of  Transverse  Shear  Stresses.  16  Layer  Hybrid,  Materials 
IV  &  V,  a/b  =  1,  hjb  =  0.3,  Simply  supported  (BC-1),  3x3  quarter  plate  model, 
1  =  0.29167,  y  =  0.29167. 
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Figure  5.1:  Eigenvalue  Coefficients  -vs-  Ply  Angle.  Case  A,  6  Layer  [+6/  —  6/  . . .],  Material 
II,  a/6  =  0.7,  h/b  —  0.1,  Simply  supported  (BC-6),  3x4  full  plate  model. 
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Figure  5.2:  Pre-Stressed  Natural  Frequency  -vs-  Ply  Angle.  Case  A,  6  Layer  [+ 9 /  —  6/.. .], 
Material  II,  a/b  =  0.7,  h/b  =  0.1,  Simply  supported  (BC-6),  3x4  full  plate  model. 
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Figure  5.3:  Eigenvalue  Coefficients  -vs-  Ply  Angle.  Case  A,  6  Layer  [+9/  —  6/  . . .],  Material 
II,  a/6  =  0.7,  h/b  =  0.1,  Clamped,  (BC-13),  3x4  full  plate  model. 
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Figure  5.4:  Pre-Stressed  Natural  Frequency  -vs-  Ply  Angle.  Case  A,  6  Layer  [+6/  —  9/  . . .] 
Material  II,  a/6  =  0.7,  h/b  =  0.1,  Clamped  (BC-13),  3x4  full  plate  model. 
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Figure  5.5:  Eigenvalue  Coefficients  -vs-  Ply  Angle.  Case  B,  6  Layer  [+0/  —  0f  . . .],  Material 
II,  a/b  =  1.0,  h/b  =  0.1,  Simply  supported  (BC-6),  4x4  full  plate  model. 
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Figure  5.6:  Pre-Stressed  Natural  Frequency  -vs-  Ply  Angle.  Case  B,  6  Layer  [+6/  —  6/  . . .] 
Material  II,  a/6  =  1.0,  h/b  =  0.1,  Simply  supported  (BC-6),  4x4  full  plate  model. 
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Figure  5.7:  Eigenvalue  Coefficients  -vs-  Ply  Angle.  Case  B,  6  Layer  [+0/  —  9/  . . .],  Material 
II,  a/6  =  1.0,  h/b  =  0.1,  Clamped,  (DC-13),  4x4  full  plate  model. 
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Figure  5.8:  Pre-Stressed  Natural  Frequency  -vs-  Ply  Angle.  Case  B,  6  Layer  [+6/  —  9/  . . .], 
Material  II,  a/6  =  1.0,  h/b  —  0.1,  Clamped  (BC-13),  4x4  full  plate  modei 
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Figure  5.9:  Eigenvalue  Coefficients  -vs-  Ply  Angle.  Case  C,  6  Layer  [+9/  —  6/  . . .],  Material 
II,  a/6  =  1.4286,  h/b  =  0.1,  Simply  supported  (BC-6),  4x3  full  plate  model. 
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Figure  5.10:  Pre-Stressed  Natural  Frequency  -vs-  Ply  Angle.  Case  C,  6  Laypr  [+9/-0/ . . .], 
Material  II,  a/6  =  1.4286,  h/b  =  0.1,  Simply  supported  (BC-6),  4x3  full  plate  model. 
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Figure  5.11:  Eigenvalue  Coefficients  -vs-  Ply  Angle.  Case  C,  6  Layer  [+ 8 /  —  9/ ...], 
Material  II,  a/6  =  1.4286,  h/b  =  0.1,  Clamped,  (BC-13),  4x3  full  plate  model. 
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Figure  5.12:  Pre-Stressed  Natural  Frequency -vs-  Ply  Angle.  Case  C,  6  Layer  [+6/ -9 
Material  II,  a/b  =  1.4286,  h/b  =  0.1,  Clamped  (BC-13),  4x3  full  plate  model 
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Figure  5.14:  Pre-Stressed  Natural  Frequency  -vs-  Ply  Angle.  Case  D,  6  Layer  \+9/ —t 
Material  II,  a/6  =  1.7,  li/b  =  0.1,  Simply  supported  (BC-6),  5x3  full  plate  model. 
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Figure  5.15:  Eigenvalue  Coefficients  -vs-  Ply  Angle.  Case  D,  6  Layer  [+0/  —  0/  ...], 
Material  II,  a/6  =  1.7,  h/b  =  0.1,  Clamped,  (BC-13),  5x3  full  plate  model. 
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Figure  5.16:  Pre-Stressed  Natural  Frequency  -vs-  Ply  A 
Material  II,  a/b  —  1.7,  h/b  =  0.1,  Clamped  (BC-13),  5 
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Figure  5.17:  Natural  Frequency  -vs-  Ply  Angle  for  All  Aspect  Ratios.  6  Layer  [+  9/— 9/  . . .], 
Material  II,  h/b  =  0.1,  Simply  Supported  (BC-6),  Various  mesh  full  plate  models. 
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Figure  5.18:  Natural  Frequency  -vs-  Ply  Angle  for  All  Aspect  Ratios.  6  Layer  [+0/—6/  . . .] 
Material  II,  h/b  =  0.1,  Clamped  (BC-13),  Various  mesh  full  plate  models. 
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Figure  5.19:  Uniaxial  Buckling  Coefficient  -vs-  Ply  Angle  for  All  Aspect  Ratios.  6  Layer 
[+6/  -  9/ . . .],  Material  II,  h/b  =  0.1,  Simply  Supported  (BC-6),  Various  mesh  full  plate 
models. 
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Figure  5.20:  Uniaxial  Buckling  Coefficient  -vs-  Ply  Angle  for  All  Aspect  Ratios.  6  Layer 
[+0/  —  0/  . . .],  Material  II,  h/b  =  0.1,  Clamped  (BC-13),  Various  mesh  full  plate  models. 
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Figure  5.21:  Biaxial  Buckling  Coefficient  -vs-  Ply  Angle  for  All  Aspect  Ratios.  6  Layer 
[+ 6 /  -  0/  . . .],  Material  II,  h/b  =  0.1,  Simply  Supported  (BC-6),  Various  mesh  full  plate 
models. 
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Figure  5.22:  Biaxial  Buckling  Coefficient  -vs-  Ply  Angle  for  All  Aspect  Ratios.  6  Layer 
[+6/  —  6/  . . .],  Material  II,  h/b  —  0.1,  Clamped  (BC-13),  Various  mesh  full  plate  models. 


Appendix  A 


Supplement  to  Equations 


A.l  Supplement  to  equations  Chapter  III 


In  the  derivation  of  the  stiffness  matrix,  eqn  (  3.64)  was  written  as 


< 


ey 


=  {<}*  =  [£]{*} 


l  J 


In  this  equation,  the  matrix  [S],  is  of  the  form: 


(A-l) 


[5]  =  [  A\  A2  A3  A4  A3  j  (A. 2) 

where 


Mi]  = 


1  0  0 
0  1  0 
0  0  1 


[a2]  = 


Tk  0  0  0 

0  ft  0  0 

0  0  r*  ft 


[a3]  = 


zak  000 

0  z/?t  0  0 

0  0  zat  z/?t 


[A4]  = 


tot  0  0  0 

0  epk  0  0 

0  0  6ak  tpk 


[M]  = 


-z  0  0 

0  -z  0 
0  0  -z 


(A- 3) 


The  column  matrix  {6}  is  of  the  form 
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mT  =  l  u«,x 


V 


o,  y 


uoyy  4"  ^o,i 


$  $  $  WOXX  WQyy  1W0'Xy  J 


where 


<Px,z  <Pyyy  *Pxyy 


Similarily,  in  eqn  (  3.65), 


-  b>  =  [SI  {«'} 


the  matrix  [St]  is  defined  as 


and 


[S,]  = 


atk  0  6' ctk  0 

o  Pk  o  6' (3k 


{6'}T  =  L  <Px  <Py 


<px 


<Py  J 


In  eqn  (  3.72): 


’1  0  0  0  0 

0  0  0  1  0 

0  110  0 
0  0  0  0  1 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  1 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  1 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

.0  0  0  0  0 


0  0  0  0  0  0’ 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  1  0  0  0 

1  0  0  0  0  0 

0  1  0  0  0  0 

0  0  0  0  0  0 

0  0  1  0  0  0 

1  0  0  0  0  0 

0  1  0  0  0  0 

0  0  0  0  0  0 

0  0  1  0  0  0 

1  0  0  0  0  0 

0  1  0  0  0  0 

0  0  0  1  0  0 

0  0  0  0  0  1 

000020. 


(A. 4) 
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In  eqn  (  3.73): 


In  eqn  (  3.76): 


'1000' 
0  10  0 
10  0  0 
0  10  0 
0  0  10 
0  0  0  1 
0  0  0  1 
0  0  10 
0  0  0  1 


[H  = 


r  [J)-1 
0 
0 
0 
0 


0  0  0  0 

[J]-1  0  0  0 

0  [7]”1  0  0 

0  0  [J]-1  0 

0  0  0  [7b]"1 


(A. 5) 


(A. 6) 


where  [J]-1  is  the  2  by  2  inverse  of  the  Jacobian  matrix,  and  [7b]-1  is  the  3  by  3  inverse 
of  the  second  order  Jacobian  matrix  (See  Appendix  B). 

In  eqn  (  3.77): 


1  0  0 
0  1  0 
0  0  [7]-1 


(A. 7) 


where  [7]  1  is  as  above 
In  eqn  (  3.82): 


[A,/t]  0  0  0  0' 

0  [a„/£]  0  0  0 

0  0  [A,/{]  0  0 

0  0  0  [a,/{]  0 

0  0  0  0  [a„/(] 

0  0  0  0  [Aw] 

0  0  0  0  [A„{j 

0  0  0  0  [A££j  . 


[■Ar,(]  - 


(A. 8) 
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where 


► 

11 

A,ij 

.  A,<  . 

(A. 9) 

[Af/r;]  —  j 

A,w 

]  (A.10) 

[A»j(]  —  | 

A.ttf 

]  (A. 11) 

[A«]  =  | 

a,€{ 

(A. 12) 

In  the  above,  the  matrix  [A]  is  as  defined  in  eqn  (  3.16),  and  the  subscripts  here  denote 
partial  differentiation  term  by  term. 

In  eqn  (  3.83): 


"  0 
0 
0 


0  (A] 
0  0 
0  0 


0 

[A] 


0 


(A. 13) 


where  [A]  is  as  defined  above. 
In  eqn  (  3.84): 


where 


and 


Here: 


[L3}  = 


I  0  0 
0/0 
0  0  B 


(A. 14) 


1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

(A  Vo) 


[B]  = 


ft  <>l 
Ol  03b 


1  0  0 
0  1  0 
0  0  1 


(A.16) 


—  ®>rjr}I'll  yiTjrjI^l 

~  V'lt  ^21 
~  ~X*Z  I'll  -  y.tf  I21 


—  X'rIV  ^12  y>r;>;^22 

H  =  ~x'rt  ri2  ~  7?  r22 

^6  =  ~x^i  ri2  -  y,Tjr;  1^22 


(A. 17) 


In  eqn  (  3.90) 
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[Sm]  =  [  Bx  B2  B3  B4  B5  ] 

where 


'  1  0  0  1  r*  0 

[Bi]  =  0  1  0  [B2]=  0  ffc 

0  0  1 J  [00 

zotk  0  da  k  0 

[B3]  =  0  zpk  [B4]  =  0  epk 

0  0  J  [00 

-2  0 

0  -2 

0  0 


(A.i8) 


(A. 19) 
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In  eqn  (  3.107): 


'1000000000000' 
0010000000000 
0000000000010 
0000100000000 
0000001000000 
0000100000000 
0000001000000 
0000100000000 
0000001000000 
0000000010000 
,  .  0000000001000 
-  0100000000000 
0001000000000 
0000000000001 
0000010000000 
0000000100000 
0000010000000 
0000000100000 
0000010000000 
0000000100000 
0000000001000 
0000000000100 


(A. 20) 


Appendix  B 


Calculation  of  Higher  Order  Derivatives 


B.l  First  Order  Derivatives 


In  relating  derivatives  in  the  local  coordinate  system  to  those  in  the  globed,  we  use  the 
standard  expressions  found  in  any  standard  finite  element  text  [146]  [9]  [23]  [136].  The 


equations  are  nothing  more  than  the  application  of  the  chain  rule.  In  finding  the  derivative 


of  a  function,  call  it  JV,  we  have: 


dN 

dN  dx  dN  dy 

(B.l) 

dr) 

dx  dy  dy  dy 

(B-2) 

dN 

dN  dx  dN  dy 
dx  df,  +  dy  d£ 

(B.3) 

dt 

or,  in  matrix  form  we  have: 


dN  > 

•  dx  dy  ■ 

r 

dr]  \ 

< 

*  ~ 

dr\  drj 

< 

' 

dN 

dx  dy 

dN 

'ST  ) 

l  3?  M  J 

{  TSy  , 

=  [J] 


dN 

dx 

SK 

Tfy 


From  this  expression  we  can  easily  get: 


dN  \ 

r  g/v 

~3Z 

1  [  =  w 

ON 

dN 

.  dy  ) 

(B.4) 


(B.5) 
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B.2  Higher  Order  Derivatives 


To  obtain  higher  order  derivatives,  we  can  apply  the  chain  rule  successively  to  the  above 
equations.  The  results  are: 


d2N  _  dN  d2x  dx  f  d2N  dx  d2N  dy\  dN  d2y  dy  f  d2N  dy  d2N  3x\ 

dy2  dx  dy2  dy  y  dx 2  dy  dxdy  dy)  +  dy  dy2  dy  y  dy2  dy  dxdy  dy  ) 

dN  d2 x  d2N  /3x\2  d2N  dxdy  dN  d2y  d2N  f  dy\2 

dx  dy2  dx2-  \dy )  dxdy  dy  dy  dy  dy2  dy2  \dy) 

and  similarity 

d2N  dN  d2x  d2N  „  32iV  dxdy  dN  d2y  d2N  /3y\2 

~d^~  ~d^W  +  !hP\dt)  +2d^d^d^  +  ~d^d^ 2  +  W\m)  (B‘7 


The  mixed  partial  turns  out  to  be: 


d2N  _  dN  d2x  d2 N  dx  dx  d2N  / dx  dy  dx  dy\  dN  d2y  d2N  dydy 
dyd£  dx  dyd^  dx2  dy  0^ dxdy  \dy  d£  +  3£  dy)  dy  dyd dy2  dy  d£ 


The  next  step  is  to  put  eqns  (  B.6)-(  B.8)  into  matrix  notation.  In  doing  so  directly  the 


result  is: 
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dN  \ 

~5* 


a3.v 

~drjr 

d2x 

dip 

(1 ?)2 

d2y 

W 

m2 ' 

d2N 

dx* 

d2N 

>  = 

d2x 

Or)  d£  \dr)  d£  or) J 

d2v 

a2N 

dxdy 

d2N 
dr)d£  J 

CD  CD 
/>>  W 
MH 

(H)2 

d2y 

w 

(1) . 

dN 

d2N 

.  dy1  , 

(B.9) 


This  equation  is  not  of  much  use  yet,  as  the  rectangular  matrix  on  the  right  hand  side 
cannot  be  inverted.  In  the  column  vector  on  the  right  hand  side,  we  know  the  first  partial 
terms  from  eqn  (  B.5)  above,  so  we  can  bring  those  to  the  left  hand  side  of  the  equation. 
This  results  in 


’  d2N  _  dN  d2x  _  d?L<Py  ' 

dr)2  dx  dr)2  dy  dr ft 

■  m2  **&  (fe y ' 

d2N 
dx 2 

d2N  d2x  dN  d2v  ON 

dr)d£  dr)d(,  dx  dr)d$  dy 

,  — 

dx_dx_  drcdriidx^dri  dy  dy 

dr)  di  dr)  d£  ’  d£  drj  orjd^ 

< 

d2  N 
dxdy 

d2N  _  dN  d2x  dN  d2v 

,  d£2  ~5x  d(,2  Hy  di2 

.  (»)2  2§?f?  (I)'  . 

d2N 

•  dy2  < 

(  d2N  'k 

Th?  I 


-  kBl 


d2N 

dxdy 


d2N 
dy 2 


(B.10) 


Solving  for  the  desired  second  derivatives,  and  substituting  in  the  terms  from  eqn  (  B.5) 
as  required  we  can  write 
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diN 

dx* 

(  d2N  02x  f-n  dN  ,  r  8N\  d*v  ( r  dN  ,  p  dN\ 

Vrila^  +  Fl2'5rJ  -  a#  [r^^  +  r^w) 

d'N 

dxdy 

Vrfi  +  3^(r2i-^  +  r22-^-) 

d2N 

-  dy1  . 

.  Sf-Si(r..^  +  r12f)-g(r2l^  +  r22^)  _ 

(B. 

Appendix  C 


BOUNDARY  CONDITIONS 

Note:  In  the  following  tables  all  boundary  conditions  are  various  versions  of  the  simply 
supported  case,  with  the  exception  of  BC-13  and  BC-14  which  are  for  the  clamped-clamped 
and  clamped-simply  supported  cases  respectively. 

C.l  Boundary  Conditions-  Predictor  Corrector  Method 

The  boundary  conditions  used  for  the  Predictor  Corrector  Method  are  given  in  Table  C.l. 

C.2  Boundary  Conditions-  Method  I 

The  boundary  conditions  used  for  Method  I  are  given  in  Table  C.2. 

C.3  Boundary  Conditions-  Method  II 

The  boundary  conditions  used  for  Method  II  are  given  in  Table  C.3. 
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Table  C.l:  Boundary  Conditions:  Predictor  Corrector  Method. 


Variable  Order:  ( ua  vQ  w0  <f>x  <py)* 


MODEL 

TYPE 

BC 

NO. 

i  =  0 

X  =  a/2 
x  =  a 

y  =  o 

y  =  b/2 
y  =  b 

1 

(01101) 

(10010) 

(10110) 

(01001) 

Qrtr 

2 

(10101) 

(01010) 

(oino) 

(10001) 

Plate 

3 

(10101) 

(10010) 

(oino) 

(01001) 

4 

(01100) 

(10010) 

(10100) 

(01001) 

5 

(01101) 

(01101) 

(10110) 

(10110) 

Full 

6 

(10101) 

(10101) 

(OHIO) 

(OHIO) 

Plate 

13 

(11111) 

(11111) 

(inn) 

(Hill) 

14 

(10101) 

(10101) 

(mil) 

(11111) 

*(1  indicates  FIXED  -  0  indicates  FREE) 


Table  C.2:  Boundary  Conditions:  Method  I. 


Variable  Order: 

(uo  Vo  lpx  <Py  W  W,x  W,y)* 


MODEL 

TYPE 

BC 

NO. 

i  =  0 

x  =  a/2 

x  =  a 

y  =  0 

y  =  6/2 
y  =  b 

1 

(0101101) 

(1010010) 

(1010110) 

(0101001) 

Qrtr 

2 

(1001101) 

(0110010) 

(0110110) 

(1001001) 

Plate 

3 

(1001101) 

(1010010) 

(0110110) 

(0101001) 

4 

(1001101) 

(1110010) 

(0110110) 

(1101001) 

5 

(0101101) 

(0101101) 

(1010110) 

(1010110) 

Full 

6 

(1001101) 

(1001101) 

(0110110) 

(0110110) 

Plate 

13 

(1111111) 

(1111111) 

(1111111) 

(1111111) 

14 

(1001101) 

(1001101) 

(1111111) 

(1111111) 

*(1  indicates  FIXED  -  0  indicates  FREE) 


Table  C.3:  Boundary  Conditions:  Method  II. 
Variable  Order:  (uQ  v0  <px  <py  w )* 


MODEL 

TYPE 

BC 

NO. 

i  =  0 

X  =  a/2 
x  =  a 

y  =  0 

y  —  b/2 
y  —  b 

1 

(01011) 

(10100) 

(10101) 

(01010) 

Qrtr 

2 

(10011) 

(01100) 

(01101) 

(10010) 

Plate 

3 

(10011) 

(10100) 

(01101) 

(01010) 

4 

(10011) 

(11100) 

(01101) 

(11010) 

5 

(01011) 

(01011) 

(10101) 

(10101) 

Full 

C 

(10011) 

(10011) 

(01101) 

(01101) 

Plate 

13 

(11111) 

(11111) 

(11111) 

(11111) 

14 

(10011) 

(10011) 

(11111) 

(11111) 

*(1  indicates  FIXED  -  0  indicates  FREE) 
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